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On the Inversion of Certain Matrices 
By Samuel Schechter 


1. Introduction. Let a; , a2, ---, dn; bi: , be,---, bp be 2n distinct, but other- 
wise arbitrary, complex numbers. For the matrix H, of order n, 


1 7 — 
(1) H ta 5} withi Si,jsn 


let G = H™* = {c;;}. (The indices i, j, k will range from 1 ton unless it is specified 
otherwise.) If, for some constant p, 


(2) a,—b;=i+j—1+7#0O, 


then H is a segment of the well known generalized Hilbert matrix, and in this 
case formulas for c;; have been given by Savage and Lukacs [4], Smith [6] and 
Collar [1]. For 


(3) a,—b;=i-—jt+p, 


Linfoot and Shepherd [3] and Collar [2] give formulas for c;; and in both cases 
Collar exhibits diagonal matrices D, K such that G = DH*K. Collar [2] and Smith 
also evaluate the quantities >> Cis bm Ci; . 

4.3 3 


These authors make use, in most cases, of the formula for the determinant [5] 
I] (a; — ar) (be — 0;) 
det H = = 
(4) “ I] (a; — bk) 
Ee 





or require the evaluation of certain involved series. 

The formulas of Collar and Smith are extended here to the general case (1). 
The method to be used does not depend on (4) but simply on Lagrange’s inter- 
polation formula. Indeed (4) comes out as a by-product of formula (17) given 
below. 


2. Formula for the inverse. Let 
(5) A(z) =[](2-a,), B(x) = [](e- 6) 
| and denote the fundamental polynomials of the Lagrangian interpolation corre- 
sponding to the a; and b; , respectively, by 


A(z) 
A’(a,)(z — a,)’ 


| where prime denotes differentiation. We then have 


B(z) 


| (6) A(z) = Boia — 6.) 


Bx) = 


Received February 13, 1959. The work for this paper was done at the AEC Computing 
and Applied Mathematics Center, Institute of Mathematical Sciences, under a contract 
with the U. S. Atomic Energy Commission. 
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THEOREM 1. The elements of G are given by 














(7) Ci; = (a; — b;)A;(b;)B,(a;) 
and if H is symmetric 
(8) Cis3 = (a; — b;)A;(b;)Ax(b;). 
Proof: For any polynomial p(x) of degree n —1 we may write 
(9) p(x) = 2) p(as) A(z) 
or 
p(x) _ C5 ‘ 
ad A@) ~ z=, 
where 
ite p(a;) 
(11) ¢ = aa,’ 
Now let p(x) = —B,(x)A(be) = p(x) and 
ae pi (ai) 
i = Aa)’ 
Then from (10) we obtain 
B,(x)A (by) oe Cri 
(12) A(z) - a-—x 


If we set x = 5; then, since B,(b;) = 4 ;, the Kronecker delta, we get that 


Chi 4 


oj = - ° 
‘ a: — b; 





Thus c;; gives the desired inverse element, that is, 


— B,(a;)A (bx) = 
A’(a:) i 
If H is symmetric: a; — b; = a; — b;, it follows that B,;(a;) = A,(b;) and 
the theorem is proved. 
We now obtain a simple formula for the sum of the c,;;. This quantity arises 
in problems of aerodynamics [2] and its simplicity allows it also to be used as a 
check on some alleged inverse. 


(a; — by) Be(a:)Ai(by) « 


i = 


CoroLiary. 
(13) cs = > (a.—bk) = 8 
tJ 
Proof. If we apply (9) to p(x) = B,(x) and set x = bh, we obtain that 
(14) > Bi(a;)Abe) = be. 


By symmetry this is also valid if A, a and B, b are interchanged. Thus we get 


that 
8 = Dia; 2) A(bi) Bias) — 2b; 2 Ai(b) Bas) = Da; — Db. 
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For the special case of the Hilbert segment (2), formula (13) gives Smith’s 
[6] formula: s = n(p + n) and for the matrix of (3) we get the formula of Collar 
[2]: s = pn. 


3. Row and column sums. Let the row and column sums of G be, respectively, 
(15) Dc = Qi, Dd cis = B; 
3 + 
and define the diagonal matrices D, , Dg by 


D. = [ay , a2, ***, On] 


Dg = [Bi , Be, >> +, Bal. 


TuEoreM 2. The matrix H then satisfies the following relations: 


ay A(b;) Ss B(a;) 
(16) aa= B’(b;) ? B; 2. A’(a;) > 
(17) H™ = D, H"D,. 


Proof: Assuming (16) to be true we note that c;; can be written in the form 


as 1 A(b;) Blaj;) a8; 
(18) ae a; BY(b;) A’(a;) a; — b: 





which gives (17). 
To prove the formula for a; in (16) we need only show that 


B(a;) =A 
™ YG b)A%e) = 


However for any function f(z) we may write 


— g(n—1) f(a;) 
(20) flas,a2, ++ ja] = 8° f(z) = ia) 
where the left side of (20) is the (n — 1)th divided difference of f(z) with respect 
to the a, (see Milne-Thomson [7] p. 9). If we apply (20) to the polynomial f(z) = 
B(x)/(x — b;) of degree n — 1 we have that 5°" f(x) = constant. Since the 
coefficient of x” in f(x) is 1, (19) is proved. 

An alternate proof of (19) may be obtained, without divided differences, by 
noting that the left hand side of (19) is the sum of the residues of the function 
B(x)/A(x)(x — b;) at the a;. (This, in fact, follows readily from (10).) How- 
ever for a sufficiently eae circle C about the origin in the complex z-plane 


B(x) 
c (x — b,)A(z) 





dx = 1 


which proves (19). 

The proof for the 8; is obtained in the same manner from (19) with the roles 
of A, a and B, b interchanged, and the theorem is proved. (The formula (19), 
incidentally, represents an extension of one of the formulas of Collar [2] for the 
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sum of a series. Extensions of other formulas given in [2] can likewise be obtained 
from (20) by specializing f(x).) 


4. Remarks. 
1.) We note that if H is symmetric then a; = 8; and if we set D = D,. = D, 
then 


(21) G = DHD 


2.) From (16) and (17) one immediately obtains the formula for the determi- 
nant of H up to a sign. That is 





(33) (det H)? = TT 
i ay B; 
and for H symmetric we get that 
(23) det H = (—1)@-n2 ll A’(a;) 
+ B(a,) 


The signs.may readily be determined, by induction, by using (6) and the formula 
for the (n + 1,n + 1) element of Gry: 


Cn+1, n4i* det Aas = det H, 


where 
. with 1 
a= b; 


/ 





IA 


H, =H and Hass ={ injsntl. 


3.) Formula (8) may also be applied to the problem of obtaining a least-squares 
fit of a function f(z) on, say, (0, 1) by a function of the form Dvr. (All vari- 


ables are here assumed to be real.) The normal equations for this problem yield 
a matrix of the form (1) with b; = —a; — 1, and the y; may be obtained by ap- 
plying G to the moments of f(x) with respect to the x. 

For fitting a function of two variables the problem to determine y;; such that 


1 1 
/ / f(a, y) — > vez ey" Pdedz = min 
0 0 t,J 
is solved by 
K = DHDFD‘H'D’ 
where 
H= — H’ = Fe Oe | K = {3} 
a; +a;+ 1)’ a! + a +1)’ ee 
if 1 1 . 
sit cenneel 
0 0 


and D, D’ are the diagonal matrices corresponding to H, H’, respectively, defined 
above. 
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Although the elements c;; get quite large in the case of the Hilbert matrix, it 
may happen that for suitable choices of the a; this may not be the case. 
An explicit solution for y; can also be obtained from (7) for the equations 


1 
i (sz) - Da") a" ae = 0, k = 1,2,---,n 
0 i 


given the moments of f(z) with respect to the at 
4.) If H is real and symmetric and if a; > b;, i = 1, 2, ---, n then it follows 
from (4) that all the principal minor determinants are positive and H is positive 
definite. In this case a; > 6; for all 7 and 7; that is, all the elements of H are posi- 
tive, since if a; < a;, then 0 < a; — b; < a; — b;. Thus B(a;) > O and, since 
A(a) has n simple zeros at the a;, A’(a;) alternate in sign. From (8) it then fol- 
lows that, if a; < a; fori < j then (—1)‘*’c;; > 0so that G has the same checker- 
board distribution of signs as the inverse of the Hilbert segment. From (16) it 
follows that in this case the a; also alternate in sign and that (—1)**"a; > 0. 
Let \ and yu denote the smallest and largest eigenvalue of H respectively, and 


let 
oe aoe (2—§). 
kk | ay = a; | 


Then it follows readily that | A.(b;) | => M7? > 1 and that 


2—2n 
» s min 2 = min ( ) 
a Ci i a; — b; 





If the a; inerease with 7 then 
min())(—1)**%e;;)* < 
. 7 
We then get for the P-condition [8] of H: 


a? 
# = max ( 





a; — i) - max ¢j; 2 max [A, (b;)]* => max Mj"~* 


so that the P-condition of H may get very large. This number has been estimated 
for the Hilbert segment (2) with p = 0 by Todd [8]. 
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A Sieve Method for Factoring Numbers 
of the Form n’ + 1 


By Daniel Shanks 


1. Introduction. Factorizations of numbers of the form n’ + 1, [1], [2], are of 
mathematical interest for at least four reasons: 1) a possible insight into the un- 
settled question as to whether there are infinitely many primes of this form; 2) a 
possible insight into the unsettled question as to whether the reducible numbers 
(explained below) have a definite density; 3) the relation of these factorizations 
to the p-adic square roots of —1; and 4) the relation of these factorizations to the 
Gaussian primes. The purpose of this paper is to describe a sieve method for factor- 
ing these numbers and to present and discuss some empirical results bearing on 
1), 2), 3), and 4) which were obtained by its use. The method can be considered 
to be based, in part, on the p-adic square roots of —1, but it is also possible to 
avoid the use of this language. 

A program based on this sieve method was written for an IBM 704 with a 
32,768-word high-speed memory, and with this program all n’ + 1 from n = 1 
to 180,000 were completely factored in about 10 minutes. Since these factoriza- 
tions of n> + 1 exceed those in existing published tables, (82 percent of these 
numbers are greater than a billion), a short summarizing statistical table should 
be of interest. In the table below, P(N) is the number of primes of the form n’ + 1 
for 1 = n S N, and for comparison, t_(N) is the number of primes of the form 
4m — 1 for1 <4m—1 SN. Further, R(N) is the number of reducible numbers 
SN,r(N) = R(N) — R(N — 10,000), and 62(N) = R(N)/N, the mean density 
of the reducibles. 


STATISTICAL TABLE OF PRIMES AND REDUCIBLES 


N r(N) R(N) 5R(N) P(N) x(N) P(N)/2-(N) 

10 ,000 2898 2898 . 28980 841 619 1.3586 

20 ,000 2935 5833 . 29165 1559 1136 1.3724 

30 ,000 2930 8763 . 29210 2268 1633 1.3889 

40 ,000 2918 11681 . 29203 2952 2117 1.3944 

50 ,000 2959 14640 . 29280 3613 2583 1.3988 

60 ,000 2912 17552 . 29253 4252 3038 1.3996 

70 ,000 2965 20517 . 29310 4888 3485 1.4026 

80 ,000 2881 23398 . 29248 5513 3933 1.4017 

90 ,000 2947 26345 . 29272 6084 4364 1.3941 

100 ,000 2875 29220 . 29220 6656 4808 1.3844 
110,000 3009 32229 . 29299 7239 5247 1.3796 
120 ,000 2934 35163 . 29303 7795 5675 1.3736 
130 ,000 2938 38101 . 29308 8369 6103 1.3713 
140 ,000 2888 40989 . 29278 8944 6531 1.3695 
150 ,000 2983 43972 . 29315 9505 6941 1.3694 
160 ,000 2952 46924 . 29328 10072 7361 1.3683 
170 ,000 2932 49856 . 29327 10658 770 1.3717 
180 ,000 2981 52837 . 29354 11223 8178 1.3723 


Received March 10, 1959. 
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SIEVE METHOD FOR FACTORING NUMBERS 79 


2. The Sieve. The sieve method (substantially more complicated than that of 
Eratosthenes) is based on the facts listed in the following theorem. Since many 
of these are well known and the rest can be easily verified, no proof need be given. 

THEOREM. For every prime p of the form 4m + 1 and every positive integer k 
there are two and only two positive solutions of: 


n< p 
4 be =0 (mod p*). 
If we call these A; and B, , then 
(2) Ar + Be = 7, 
and if 
(3) h = }3(p+1)A;i (mod p), 
and 
(4) C. = h(Ay’ + 1)/p" (mod p), 
the A; for k = 2, 3, --- may be computed recursively from A; by 
(5) Aras = Ax + Crp. 


Further, for all positive n, 
n+1=0 (modp’), 
if and only if n is given by one of the linear forms: 
A, + mp* 
(6) ders k (m = 0, 1, 2, ---), 
B, + mp 


and aside from these factors of n* + 1 (obtained as p runs through all primes of 
the form 4m + 1) the only other prime-power factors are the obvious 


(7) n’+1=0 (mod2) forn =1 (mod 2). 


We will adopt the convention that A, is the smaller of the two roots of (1) 
for k = 1, and note that this implies: 


(8) Ai < 3p < B,, 


but it does not imply that A, < B, if k > 1. For example, let p = 5. Then A; = 2. 
Thus h = 1, and for k = 1, 2, 3, --- we compute: 


Ax = 2, 7, 57, 182, 2057, 14557, --- 
B, = 3, 18, 68, 443, 1068, 1068, --- 


lI 


and note, As > B;. 

We also have in this case, the interesting degeneracy B; = By, . While similarly, 
for p = 13, we have the degeneracy A; = A,, it can be seen from (1), and it is 
important for the validity of the sieve method, that for every p: 
™ < Ay 


B, < B:. 


(9) 
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Now suppose we wish to factor every n° + 1 for 1 < n < L. The sieve method 
proceeds as follows: 

A.) In L contiguous cells we write the corresponding values of n’? + 1. 

B.) We now divide out the factor of 2 for each n = 2m + 1, (m = 1, 2, ---) 
and store the quotients back into the same cells. 

C.) At n = 2, we find the number 5. This implies (see a general proof below) 
that 5 is a p and, corresponding to it, A, = 2. Using formulas (2) through (6), 
we factor a 5 from each n’th cell, where 


(2 + md (m = 1, 2, ---) 
|3 + ms (m = 0, 1,2 


r= 


and store back the quotients. Then a second 5 is factored for each 
[ 7 + m25 
|18 + m25 


and so on until both A; and B, are greater than L. By now all factors of 5 have 
been removed. 

D.) At n = 3, we now find the quotient 1. This means 3 is not an A; number. 
It is, in fact, reducible (see below). We so record it and move on. 

E.) At n = 4, we find 17 and divide out all factors of 17 as in C. 

F.) At n = 5, we find 13 and divide out all factors of 13. 

G.) Proceeding in this way we examine the contents of the n’th cell when we 
get to it, and find one of two cases. 

1. The quotient is 1 and no new prime factor is contained in n’ + 1. 

2. The quotient is >1. In this case the quotient must be a new p = 4m + 1 
and n must be its corresponding A; number. Proof: Since the quotient cannot 
have a factor of 2 or a prime of the form 4m — 1, every prime factor must be a 
p = 4m + 1. But for such a p, n must be A, for, if p had occurred in an earlier 
n’ + 1, its A, would have been encountered earlier and p would have been fac- 
tored out. By (9), n’ + 1 is not divisible by p’. Nor can n be simultaneously an 
A, for two distinct primes, since from (8) the quotient: 


(A, + 1)/p S$ 3A, 


and therefore cannot contain a second prime which (again by (8)) must be >2A),. 

This completes the proof, and shows we are obtaining complete factorizations 
without the use of any trial divisions and without the need to know in advance 
the primes of the form 4m + 1. They are, in fact, being generated (though not 
in numerical order) by the process itself. 


r= (m = 0, 1, 2, ---) 


3. The Program. In the 704 program mentioned above the following changes 
were desirable or expedient. 

a.) The sieve was started with —(n’ + 1) in the n’th cell and the absolute 
value of each quotient was stored back in. If no division occurred in a particular 
cell the contents of this cell remained negative and this indicated that n> + 1 was 
a prime. 
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b.) Since in a binary machine division by 2 is so simple, step B was combined 
with step A. 

c.) The factoring out of any prime p began when n reached its B, instead of 
its A; . This saves time and memory and is equally valid. 

d.) The linear forms, (6), are easily obtained by use of an index register. The 
high efficiency of this method makes possible over 3000 divisions/sec. 

e.) The sieve is done in several segments. In a memory of 32,768 words, with 
4096 words reserved for an operating routine, a sieve segment of 20,000 words is 
possible, since nearly 7000 words (see g) must be reserved to save the primes (and 
corresponding A, numbers) for later phases. The program itself, although quite 
intricate, takes less than 300 words. 

f.) An upper limit, L, of 180,000 (9 segments) was chosen since a tenth segment 
would cause n® + 1 to overflow the length of a 704 word. This is , (nearly 34.36 
billion). Of course n could be made larger by using double precision. 

g.) Although 127,162 of the 180,000 numbers are A; numbers, in most cases 
the corresponding p is enormous and is not used in sieving—the criterion being 
p + A, S L. The number of p’s actually used was 6693. 

h.) The program was checked by being repeated with a segment of 20,500 
This changed all the phase relations in the linear forms (6). It also raised L to 
184,500 and we find P(184,500) = 11,486; R(184,500) = 54,162. 

i.) If it were not for the “degeneracies” mentioned just before (9), a much 
simpler sieve method would be possible. One would not compute the A, and B, 
sequences and take out all factors of p at once but could instead re-factor the 
Aog+ mp’, the B, + mp’, the A; + mp’, etc., when n came to Az, B., A;, etc., 
respectively. This was in fact attempted, and some very rare errors ensued. For 
example, one found R(20,000) = 5832 instead of the true 5833, the single error 
stemming from the degeneracy B; = B, (for p = 5) mentioned above. 

j.) The main output of the program as described was the table shown above 
but in a much greater detail (the interval in N being 100 instead of 10,000). The 
program actually computed 7,(N) by counting the p’s instead of the x_(N) shown. 
A small modification of the program also produced a 360-page table showing the 
largest prime factor in every n° + 1 for n = 1 (1) 180,000. Small auxiliary tables 
concerning the p-adic square roots (the A, and B, sequences above) and the dis- 
tribution of the reducible numbers were also produced. 


4. The Primes. Now consider P(N) in the statistical table. We see a steady 
growth which is nearly proportional to r_(N ). Since 


N 
1r_( N ) ww 1 de 


2J2 Inz 


by the prime number theorem, the numbers P(N) are in excellent agreement with 
the conjectured formula [3] of Hardy and Littlewood: 


N 
(10) P(N) ~ 0.68641 [ 
2 nz 


which therefore implies P(N )/x_(N) ~ 1.3728. 
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The constant in (10) is equal to the infinite product, taken over all odd primes, 


sAL-G)e-»], 


where (—1/p) is the Legendre symbol. It was computed by A. E. Western [4] 
who also verified the substantial correctness of (10) to N = 15,000. Since the 
infinite product converges too slowly, Western used a transformation of the product 
due to Littlewood. The following related formula is simpler: 


v= 3 $6) ( 5247) ( ~ ap) 
(10a) 0.68641 4 Gt) ate 1+5> 1 op — 1%): 


Here G is Catalan’s constant, and the product is taken over the 7, primes. More 
rapid convergence may be obtained by multiplying through by the identity: 


_ [17 _5(8)_ eoaht 
(10b) ee E r@r@ >t, (: +e ”) 


where L(4) = ia (—1)” (2n + 1)“. With this improvement, the first two p’s, 
i.e., 5 and 13, already suffice to yield the five decimal places shown. 

Among the Gaussian integers, a + bi, the Gaussian primes on the positive real 
axis are those of the form 4m — 1, and thus are those counted by z_(N). On the 
other hand, n + 7 is a Gaussian prime if and only if n* + 1 is a rational prime. 
Therefore, the column headed “P(N)/z_(N)” shows that the +1 horizontal 
line in the Gauss plane (and therefore also the —1) has about a 37 percent greater 
density of primes than the axis has. (For an attractive picture of these primes see 
van der Pol’s Tea Cloth, [5].) 

It is of interest to mention briefly the ‘‘Gaussian twin’ primes, i.e., those where 
(n — 1)? + 1 and (n + 1)” + 1 are both prime. They are not at all rare. In the 
last block of 1000 numbers with L = 184,500, we find no less than 9 pairs, namely: 
n = 183585; 183635; 183685; 184055; 184075; 184145; 184185; 184325; and 184495. 
The last pair, 34038036037 and 34038774017, were the largest primes obtained 
by this program. The conjecture suggests itself that there are infinitely many 
such “twins.” 

It may also be mentioned that (although they were not particularly sought) 
the program yielded a very large collection of ‘large’ primes, i.e., those over 10 
digits. Not only are 4830 of the n’ + 1 “large” primes, but many others are equal 
to twice a “large” prime. For example, 184499’ + 1 = 2-17019940501. 


5. The Reducible and Irreducible Numbers. A reducible number, r, is a positive 
integer whose arctangent is a linear combination, with integer coefficients, of the 
arctangents of smaller positive integers: 


r—1 
(11) tan’ r= >a, tan’ n. 

n=l 
If no such linear combination is possible, the number is irreducible [6]. For example, 
since 


tan’ 3 = 3 tan” 1 — tan’ 2, 
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3 is reducible. So are 7 and 8, while 1, 2, 4, 5, 6, 9, and 10 are irreducible. (If we 
replace the arctangent function in (11) with the logarithm, we would be defining 
composites and primes instead, so that the irreducibles can be thought of as the 
quadratic analogue of the primes.) John Todd showed [6] that r is reducible if and 
only if the greatest prime factor of r’ + 1 is less than 2r. For each irreducible, n, 
(other than 1), there is thus a unique prime p > 2n which is the greatest prime 
factor of n’ + 1. [6, p. 526]; and this is the 1 — 1 correspondence indicated above 
between A; and p. A; numbers are therefore irreducible, and all others except 1 
are reducible. This includes all Az, A;, ---, ail B,, Bs, ---, and some numbers, 
such as 21, which are neither A, nor B; . 


6. The Density of Reducibles (Heuristic). If R(N) is the number of reducibles 
<N, and 6g(N) = R(N)/N, and if Limy.« 5g(N) exists, we call this limit the 
density of the reducible numbers, and write it dz . 


(12) de = Limy.« R(N)/N, if it exists. 


However, its existence is still unsettled, as is the older, and somewhat related 
question: 


(13) Does P(N)—~ ? 


Early counts of the reducible numbers up to 5000, by J. C. P. Miller and Todd, 
[6], [7], [8], were based on Todd’s criterion [6] and Wrench’s table [2], and indicated 
that 52(N) persisted in the vicinity of 0.29. The present paper extends N to 184,500 
and shows a continuance of this state of affairs with a slow growth (in the mean) 
to about 0.293. 

Chowla and Todd [7] have proved that if C(N) is the number of numbers 
n < N for which the greatest prime factor of n is greater than 24/n, then the 
Limy..C(N)/N exists and equals In 2 = 0.693---. But note that here n ranges 
through all numbers, not just those of the form m* + 1. Jf this latter class of 
numbers constituted a good sample of all the numbers as regards factorization 
properties, (with a rigorous definition!) it would follow from 24/m? + 1 ~ 2m 
and Todd’s criterion of reducibility that 5g does exist and is equal to 1 — In2 = 
0.30685: - -. 

But, if such a deduction were possible, it would probably also follow from the 
same line of reasoning (good sample concept) that the mean local density of primes 
of the form m’ + 1, like that of the ordinary prime sequence, would be 


1 1 
iIn(m? +1) 2inm 





and thus 
‘ 1" dm 
PN) ~3/, In m w_(N). 


Now in fact we have seen that P(N) remains persistently and significantly higher. 
And this greater prevalence of primes is consistent with the smaller fraction of 
reducibles (.293 instead of .307) which is observed. It is, of course, not precluded 
that 52(N) will rise to .307 for very much larger N. Its slow growth, mentioned 
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above, seems to be of order of a/In N and to be associated with the falling off in 
the mean local density of the excess primes, P(N) — 2_(N). This (average) 
increase of 52(N) hardly shows in the Table since In N changes so slowly and the 
fluctuations are almost of the same size. Nonetheless it is there and may carry 
52(N) nearly up to 0.307. 

The gross facts on which any attempted assessment of the “‘good sample” 
concept must be based are these. The n” + 1 numbers have only about } of all 
primes, [9], as possible factors, but in “‘compensation” these factors occur twice 
as often (see (6) above). While this may suggest a factorability of the same order 
of magnitude, certainly no more exact equivalence is implied. 

Assuming the future establishment of R(N) ~ 6z-N, the (deeper?) question 
of O(R(N) — ée-N) will arise. Concerning this question—that of the uniformity 
of the distribution of the reducibles—the following table and figure are informative. 
Each of the 1800 intervals of 100 numbers: 


100m <n < 100(m + 1) (m = 0,1, ---, 1799) 


has at least 17 reducible numbers and at most 42. The 52,837 reducibles < 180,000 
are distributed as follows: 

r,no. of reducibles 17| 18| 19| 20|21|22|23|24|25| 26| 27| 28| 29 

v, no. of intervals 5| 6| 5| 13|22|29|69|82| 93] 127| 154|147|167 





30| 31| 32| 33|34|35|36|37|38| 39] 40| 41| 42 
179|140|147|130|93|72|41|31|16| 16] 11] 1| 4 





In the Figure a bar graph of this distribution is compared with a binomial dis- 
tribution, »v (r): 

- 100! r ; 100—r 
where the “probability,” 0.2935, was taken to be the final value of the mean 
density, 52( 180,000). 


7. P-adic Numbers and Degeneracy. The sequence of the partial sums of the 
infinite series: 


(15) A,+ 2 CoP’, 


where p is a prime of the form 4m + 1 and A, and the C; are determined by (1) 
through (4), is a convergent sequence in the p-adic sense. (See [10] for an elementary 
account of p-adic valuation.) Since it converges, it represents a p-adic number, 
namely, one of the two values of ~/ —1. This sequence, it is seen from (5), is the 
sequence A; . The other sequence, B, , converges in the p-adic sense to the other 
~/ —1 and their sum, A; + B; , converges in the p-adic sense to 0. 

As an example, take p = 5, and A; = 2, 7, 57, 182, etc., as above. If we write 
them in the quinary system we see that they represent the sequence obtained 
from the 5-adic number 


(16) - + +3140223032431212. = ~/—1 
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The Continuous Curve is the 
Binomia! Distribution: 
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Fig. 1.—Distribution of the reducible numbers between 1 and 180,000 into the 1800 in- 
tervals of 100. 


by starting at the quinary point and taking more and more places to the left. If 
we take k places and carry the computation to at least k places we find that 
A? + 1 = ----000---00. 


k zeroes 
so that A,’ + 1 is divisible by p* and A, is an approximation to ~/—1 correct 
(p-adic sense) to k places. Similarly the B, sequence gives the complement: 


(17) - + + 1304221412013233. = — /—1. 


The point of this review is to indicate, first, that the degeneracy B; = By = 
1068 mentioned above (9) is associated with the last zero digit in (17). That is, 
B; = 13233 (quinary) and so is B,. Similarly from the zero digits of (16) we find 
that A, = As and A, = Ay. Since an irrational p-adic cannot have periodic 
digits, [10, p. 196], it is reasonable to conjecture that the digits are “normal.” 
This would imply that degenerates appear infinitely often in every A, and B, se- 
quence for every p, and further that we should expect multiple degenerates, A, = 
Ars: = Azye, ete. For all that, aside from the B; = B, = 1068 for p = 5 and 
the A; = A, = 239 for p =.13 already mentioned, there are no other degener- 
ates = 180,000. 


8. The Difficulty of the Unsettled Questions. It has often been remarked that 
questions like (13) are “very difficult.”” The intent of this section is to assess this 
difficulty. We do this by comparing the very simple Eratosthenes sieve for the 
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ordinary prime sequence with the present one for n’ + 1 and note that the latter 
is more complicated in the following three ways: 

1.) Instead of one linear form, mp, for each prime, we have a double infinity: 
Ay + mp*, Be + mp". 

2.) Instead of a zero origin we have A; and B, origins, which are not related to 
p in any simple fashion. While A, , A;, --- and B,, B,, --- can be computed by 
the more or less complicated relations (2) through (6), A: can arise at any n 
satisfying ~/p — 1 S n S (p — 1)/2. Further, we have the complications of 
occasional degeneracy. 

3.) Finally, whereas in the Eratosthenes sieve it is not necessary to divide, but 
it suffices to scratch the cells in the linear form, here we must divide the p out. 
Otherwise, we would not obtain the new prime hidden in each A,’ + 1 which is 
not itself prime. 


9. Generalization. In conclusion, it should be stated that while we have con- 
fined ourselves here to n’ + 1 the same type of sieve method is applicable to 
n’ + a, for a = +2, +3, etc. The main change is that these programs would be 
based on the p-adic square roots of —a The sieve itself would generate those 
primes for which —a is a quadratic residue, that is, all possible divisors. It is thought 
that comparable statistics will be found for the primes of these forms and for the 
“generalized irreducibles,” that is, those n which yield a new prime factor. This is 
now being investigated, [11]. 
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Numerical Quadrature of Fourier Transform 
Integrals II 


By H. Hurwitz, Jr., R. A. Pfeiffer and P. F. Zweifel* 


1. Introduction. In a previous paper [1], even-point Gaussian integration 
formulae for the numerical quadrature of the trigonometric integrals 


(1a) S(2z) =| $(k) sin kx dk 
0 

and 

(1b) C(x) = f: W(k) cos ke dk 
0 


were presented. In this note we present odd-point Gaussian formulae. In addition, 
by making use of the fact that in both the even and odd point schemes the points 
are equally spaced, we derive a method for evaluating the integrals in Eq. (1) 
for a large number of values of the parameter x from a knowledge of the functions 
¢(k) or ¥(k) at a specified set of points {k,}. In other words, we remove, to some 
extent, one of the disadvantages of the method described in [1], to wit, that for 
each different value of x the functions y and ¢ had to be evaluated at an entirely 
new set {k,}. 

This disadvantage is also inherent in a scheme proposed by Goldberg and Varga 
[2], whose method also appears to be somewhat more difficult to apply than that 
used here. 


2. Odd Point Gaussian Formula. In [1] the integral (la) was transformed into 
the following form (the arguments for (1b) are identical) : 


4 
(2) S(z) = = [ cos ryo(xz, y) dy 
2x Ly 
with 
‘ = 2 1 
(3) o(y,z) = Do (-y'6(2[v +44) 
The (2N point) Gaussian formula gives for the integral 
N (») . 
(4) S(z) = 2 2M acy”, 2) 
2 j=1 cos ry} 
where the y;“” are roots of the equation 
(5) Tews (cos wy) _ 9, 
cos Ty 





Received June 12, 1958. 
* Present affiliation: Department of Nuclear Engineering, University of Michigan, Ann 
Arbor, Michigan. 


87 











88 H. HURWITZ, JR., R. A. PFEIFFER AND P. F. ZWEIFEL 


Here T2y +1 is a Chebyshev polynomial of the first kind [3]. The W; are the Christof- 
fel numbers which may be obtained from Eq. (14) of [1] which was derived from 
the fact that the Gaussian formula gives a rigorous answer if o(y, x)/cos ry is any 
polynomial of degree up to and including N — 1 in cos’ zy. 

The 2N + 1 point formula N = 1, 2,3 ---- is similar to Eq. (4), except that 
the y;°” are the roots of the equation 


(6) Uon41 (sin wy) = 
where U,, is a Chebyshev polynomial of the second kind [3], 
sin (n + 1)xy 
(7) we (cos ry) = = — sney 
Then, the quadrature formula is simply 
N (N) 
(8) S(e) = £) W8o(0, 2) + 2 aly, x)} 
j=1 CO wy” 


where the first term of (8) has been split off since Eq. (6) has a root at y = 0 for 
all N. 

The Christoffel numbers are found from the following (N + 1) simultaneous 
equations: 


1 TA +3) _ f ow —2 Nr «| 
(9) wrath | +2 cos ” Yj W; 7 
The y; (for the 2N + 1 point formula) are 
(10) yo” = j= 1,2,3,--°°N. 


Fig. 1 of reference [1] showing the error in percent in the 1, 2, and 4 point formulae, 
is reproduced as Fig. 1 of this paper, with additional curves showing the errors in 
the 3 and 5 point formulae. The summation procedures are, of course, identical 
for odd- and even-point formulae. 


3. Integration Over a Fixed Point Set {k,;}. In many applications, it may be 
necessary to evaluate the integrals of Eqs. (1) for a large number of values of the 
parameter x. In general, for each x, the transform function must be evaluated at 
a different set {k;}. However, by taking advantage of the fact that the points of 
both the odd and even point Gaussian quadrature schemes are equally spaced, it 
is possible to choose initially a set {k,;} at which the function ¢(k) is to be evaluated, 
and these points are used in the integration scheme for each value of x. Clearly, 
S(x) cannot be obtained for all x; one initially chooses a value %max which is the 
maximum value of the argument for which one wishes to evaluate the function 
S(x). Then it is possible to find S(2max/}m), from the ¢(k;), m = 2, 3, 4 ----,* 
which in practice, will frequently be enough points to allow a smooth curve of the 
functions S(x) to be drawn. 


* For the cosine integral, m = 1 is also a permissible value. Thus, tmax may be chosen as 
one-half the maximum desired value of zx for which the cosine integral is desired. 
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Fig. 1.—The per cent error introduced by evaluating the integral C(z) by the one, two, 
three, four and five point formulas, with ¥(k) = (1 + L*k?)~ plotted as a function of z*; z* = 
2/L. 
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additional points are used for evaluation of the smaller values of x. For each such 
value, some of the points will be “wasted”, since the k; are chosen such that the 
zeros of sin(kx) always fall at one of the k; for the values of x that we consider. 
We make this choice because it appears to be the most convenient for our purposes. 
It means that for every n, S(2max/(n/2)) can be evaluated with equi-spaced points, 
the points falling at the zeros of sin ktmax/(n/2) being weighted with a Christoffel 
number W; = 0 

Fig. 2 illustrates the points {k,} at which ¢(k) is evaluated, and the possible 
values of x which may be found from these {k;}. 

It is clear, from the figure, that the integral involving sin kamax can be evalu- 
ated from a 1 point formula, that involving (sin ktmax)/(3/2) from a 2 point 
formula, that involving (sin krmax)/2 from a 3 point or a one point formula, etc. 
If convenient, for machine calculational purposes, the end points can be included, 
so that an N point formula becomes an N + 2 point formula with the Christoffel 
numbers of the first and last points of each half cycle being set equal to zero. This, 
as was mentioned previously, makes the point spacing equidistant everywhere. 

Consider S(2max/(n/2)). Clearly, not every N point Gaussian formula involving 
the points {k,;} can be used. We find the available formulas in the following manner: 

We write all the factors of n including n but not 1 (not only prime factors). 

Then, if these factors are fi, fe, f;, the available formulas are (f; — 1) 
point, (fe — 1) point, (f; — 1) point. (Note that when n is prime only one and 
(n — 1) point formulae may be used.) 

For an (f; — 1) point formula, ¢(k) is evaluated at the points k, , where 


a j = 0,1, 2, ----, fi); (fi, fe + 1,--++ 2fs); (2f:, 2fs + 1, ---+ 3fi). 
‘ first half cycle second half cycle third half cycle 


The Christoffel numbers associated with the points corrésponding to 7 = 2, ---- 
fi—1;f: +1, ---+ 2f; — 1; ete. are those obtained for the (f; — 1) point formula, 
while the end points hf; , h = 0, 1, ---- have Christoffel numbers zero in order to 
simplify machine coding. 


General Electric Research Laboratory, 
Schenectady, New York and 

Knolls Atomic Power Laboratory, 
Schenectady, New York 


1. H. Hurwitz, Jr. & P. F. Zwetrex, ‘Numerical quadrature of Fourier transform inte- 
grals, ™ - 55, MTAC, v.22, =v. 140-149. 
R. Gotppere & R ARGA, “Moebius inversion of Fourier transforms,’’ AEC 
Report WAPD. TN-524, 1956. 
A. Erpe yi, W. Maanvs, F. OperRHeEeTTINGER & F. G. Tricomi, Higher Transcendental 
pa, v. 2, McGraw Hill, Inc., New York, 1953. 


et 


Dr. 


~~ = 











~~ ra. ~~ aa | 


EF Fei de 





LN vn = AEE Bie Sd 





Note on Bivariate Linear Interpolation for 
Analytic Functions 


By Walter Gautschi 


1. Introduction. As more tables of complex analytic functions are being pro- 
duced* the problem of complex interpolation has recently obtained increasing at- 
tention. The research has mainly been directed toward formulas of high accuracy 
based on Lagrangian and Lagrange-Hermitian interpolation (see, e.g., [9], [10], [11] 
and the literature quoted there). In the present note we are dealing with the simpler 
means of bivariate interpolation, having an occasional computer in mind with only 
moderate accuracy requirements. It will be shown that bivariate interpolation 
enjoys special properties if applied to the real and imaginary part of an analytic 
function and we shall indicate possibilities for slightly improving the accuracy. 
Although we are restricting ourselves to linear interpolation, similar arguments 
could be used for bivariate quadratic interpolation. 


2. Bivariate cartesian interpolation. Let f(z), z = x + iy, be a single-valued 
analytic function in a domain D of the complex plane. We consider first bivariate 
interpolation with respect to a rectangle R in D. 

By a translation and rotation of the z-plane we can fix R to be the rectangle 
|x| S $h,|y| S 3k for some h > 0, k > O with h = k. We shall adopt the nota- 
tions 


(1) 2p.q = >ph a 3tgk, So.0 ons S (29.4); 


so that points z,,, in the closed rectangle are characterized by |p| S 1,|q| < 1. 
We also set k = ch(o S 1). 

Given the values of f(z) at the vertices of R, bivariate linear interpolation to 
evaluate f,,, proceeds as follows. First we interpolate linearly along two parallel 
sides of R, say the sides parallel to the z-axis, to obtain 


fpa=41-pifaanthlt+rp ha, fea = 41 — p)ifaat 1 +t p)fa- 
Then we interpolate linearly in the y-direction between these two values: 


fre = 311 — Qfpathlt+ Ofsa 
Thus 


(2) foe=HO- pl —- fant (14+ pl - Of 
+ (1-—p)\Qlt+q@)faat (14+ p)(1 4+ gif,} 


In order to investigate the error fp,¢ — f>,¢ We expand this difference into powers 
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of h by means of Taylor’s series. We find 


Hepa) — fou = & (PP - Ff — 1+ 0°)f"(0) 


‘ 3 
(3) +h {p — 30 pq + (30° — 1)p + iol3p'q — og — (3 — o ql }f’”’(0) 
+ O(h*). 
From (3) we note the following facts: 

(i) Bivariate linear interpolation based on the vertices of a rectangle R with 
side ratio o has an error of order h” except for points zp,, lying on the 
hyperbolas p’ — oq — 1 + o° = O where the error is of order h’; 

(ii) the interpolation error is largest in absolute value in the neighborhood of 
the midpoints of the larger sides of R; 

(iii) in case of a square (o = 1) the interpolation error is of order h’ except for 
points along the diagonals, where it is of order h® and the center of the 
square, where it is of order h’. 
At the center of a rectangle the result of bivariate linear interpolation is simply 
given by the arithmetic mean 


(4) fas = (faa + fia t+ faa + fia). 


The relation (4) in case of a square is often used for the numerical solution of the 
Laplace equation. In this connection the error property mentioned last in (iii) 
is a well known fact. 

The magnitude of the coefficient of h’f’(0)/8 in (3) is shown in Figs. 1 and 2 
for the cases of a rectangle (¢ = 4) and a square (o = 1). 

3. Successive bivariate interpolations. Suppose now that f(z) be tabulated at 
the nodes of a square grid with grid width h. The simple device of forming the 
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Fig. 1.—Curves of constant interpolation error in a rectangle (neglecting third order 
terms). 
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- Fic. 2.—Curves of constant interpolation error in a square (neglecting third order 
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Fig. 3a. Grid line. Fic. 3b. Grid line. 
1 : arithmetic mean (4) yields the value of f at the center of a grid square accurately 
to 0(h*). Therefore we can generate, within this accuracy, a new table of f with 
i respect to the ‘“‘dual” grid consisting of all centers of the original grid. With the 
union of both grids as a new starting grid the process can be repeated indefinitely. 
Thus, in principle, it is possible to evaluate f(z) for any z with an error of essentially 
0(h*) by a succession of linear bivariate interpolations. 
For the purpose of interpolation, of course, one should try to introduce new 
; grid points as economically as possible and in such a way that the given point z 
- , will be close to the center of a new grid square. 


To illustrate, assume that z lies near the midsection of a grid line ab (see Fig. 
3a. ). 
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In this case we evaluate f(m,), f(mz) at the centers m, and m, by averaging and 
then apply bivariate linear interpolation to the square am,bm, . It should be noted 
that the parts with respect to this square are 


la — ta | |a — us| 
, SO = _ , Se ————  —-SC—i‘— 
Pr la — ™ | Pi P2 ; P2 [a — ma | Pi + pr, 
if p1 , P2 are the corresponding parts with respect to the square abcd and p; > p. > 0. 
In a case like that shown in Fig. 3b we may first compute f(m), f(mez), f(ms) 
by averaging and f(m), f(m2) again by averaging. Then we can use the square 
mbn2m for bivariate interpolation. 


4. Bivariate polar interpolation. Let f(w),w = pe”, be a single-valued analytic 
function defined in a certain domain EF of the w-plane. We consider now bivariate 
interpolation with respect to a sector S of an annulus located in EZ. Without 
restricting generality we can assume that S is the sectorr S$ p S r + Ar,|@| < 
3A0@. In analogy with (1) we denote 

p+1 


(5) trp = (r+ PE? ar) exp (Hig 8) (lp| S1,/¢| 81), 


(6) fo.a - F(Wp,¢); 


and also set A@ = rAr. 

With these notations the formula for bivariate polar interpolation is identical 
with formula (2). Its properties are readily deduced from the results of section 2 
by means of the conformal mapping z = In w. This function maps S onto the 
rectangle R{ln r < x S In(r + Ar), | y| S $46}, the point m = +/r(r + Ar) of 
the w-plane being carried over into the center of R. The image of w,,, is given by 


(7) 254 = (Inm) + 3p~h+ 3iqk = * 


where 





(8) 
2, Ar 1 Ar 
= p+ (1p) {1-2 (8-+29) 
1 2, (Ary \ 
+3(7 + 8 +39") () +---}, 
(9) q=4% 
(10) h=in(1 +), k = AQ = rAr. 


* For the sake of simplicity we use here the same notation as the one introduced in (1) 
which refers to the special case m = 1. 
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Now let 
f(w) = fle) = F(z), F(zas) = Fas, Fes = }{(1— @)(1 — 8)Fa 
+ (1+ @)(1 — B)Fi4 + (1 — a)(1 + 6)P in + (1 + @)(1 + 6)Figl. 
Then it follows from (3) 
F(254) — Fea = (ph — 7k —W + 2)F" (Inm) + --- 
or, in view of (8)—(10), 


(11) Flzs) - F< = 3(*) (p — rd —1+779)F’ (Inm) + O[(4r)4. 


Since furthermore 


Fr. — Fou = fa? {-Q-—@Faat(-@Fia- (0+ @PFaat 1+ Fi} 





—_ 2 3 
wt 5 P (*) F’ (In m) + O[(Ar)’, 
we have 

F(25,4) — Fy. = (F355 — Fo.) + (F(25.a) — F3;) 


- (7) {(1 — p)F’ (In m) + (p — o'r’ — 1+ 7°P°)F” (In m)} + Of(r)’). 


Thus, if we transform back to the w-plane, observing that 


F’ (In m) = mf'(m), -F” (In m) = mf'(m) + m'f"(m), =m = r+ 0(A4r), 


we get 
f(wp.e) — So.0 
_ (ar)? POD ae a , 2 222 22\ en \ 3 
“fo q )rf'(m) + (p — rrg —1+78r)f’(m)} + Of(ar)’). 


In order that the coefficient of (Ar)* be zero for all choices of f(z) we must have 
|p| = |q| = 1. Hence, in general, linear bivariate interpolation for all points of S 
is affected with an error of O[(Ar)’] and there are no exceptional points other than 
the vertices of S with higher order accuracy. 

We note, however, that exceptional points will be present if we adopt a modified 
version of bivariate interpolation, namely, using f as defined in (8) in place of p. 
It follows then from (11) that the interpolation error becomes of order (Ar)* for 
points w,,, along the curves 


p —tr¢d —1+7rr =0. 


Moreover, if A@ = In (1 + (Ar/r)), that is, if the conformal image R of S is a 
square, this modification leads to an error of order (Ar)* at the point 
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m= +~/r(r + Ar) of S. These facts may occasionally be used to improve the 
accuracy of bivariate polar interpolation. 


American University 
Washington, District of Columbia 
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12. Added in Fah In a private communication E. Kreyszig and J. Todd have studied 
the case of bivariate quadratic interpolation with respect to nine pivotal points zp,.¢ = ph + 
igk (k = xh and p, q running independently through —1, 0,1). For the interpolation error 
they obtain the expression 


P hs 
f (ena) — Soe = [ vlot — 1) — tg (@# - 1) |F 6"" + WX" - 1) 
+ (qt — 1) + Sixpgtp — 1 — (gt — NE fer + 010), 


from which they observe that no exceptional points exist, other than the pivitol points, 
where the error is of order A‘. 
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Gaussian Quadrature of Some Integrals 
Involving Airy Functions 


By J. A. Riley and C. Billings 


1. Summary. The integrals 





, Oe tan | , AB i .B 
a [ Fa pe ; = pt : x p , s pe 
one 4 [ , A” Te f , B 
[ zx pe A zx F2 dx, a Hr rr dx, q zx pat 


(in which A, B are the Airy functions Ai, Bi, and in which A’ = dA/dz, B’ = dB/dz, 
F’ = A’? + B’, F” = A” + B”) are the coefficients in certain power series arising 
in the diffraction theory of electromagnetic waves. The values of these integrals 
i were computed for \ = 0, 1, --- , 20 by dividing the interval of integration into 
four closed subintervals (0, 8), (8, 16), (16, 24), (24, ©) and carrying out the 
integration over each interval separately. The integrations in each of the first three 
’ } intervals were performed by using a 27-abscissa Gaussian quadrature scheme. 
The integration over (24, ©) was effected by exact integration of an asymptotic 
form for the integrands. The values of the integrals are shown in Table 1. Although 
an exact estimate of the total error was not obtained, it seems reasonable that the 
results are correct in most cases to five significant figures. An independent check 
r on the accuracy has been provided by certain calculations made at the Electro- 
} magnetic Radiation Laboratory, Air Force Cambridge Research Center, in which 
our results were used. 





2. The Integration Scheme. The integrals (1) have the general form 
(2) (A) = [ af(x) da. 
; 0 


) An examination of the values of the integrands will show that they have a skew 
‘bell-shape’ rising from zero at x = 0 to a maximum point somewhere in the in- 
terval (0, 8), and thereafter decaying quite rapidly to zero. On the basis of this 
behavior, and for convenience, the interval of integration is subdivided into four 
intervals: (0, 8), (8, 16), (16, 24), (24, ©). In the course of the numerical evalu- 
ation, it is easy to check whether sufficient accuracy is obtained by retaining only 
the first integral, or the first two, etc. Thus we write: 


STI 5 elena te 


hans 


(3) IA) = [ a f(x) dx +f a f(x) dx + [. af(x) dx + [ a f(x) dz. 


6 a a ie 


Received August 14, 1958. The work reported here was done under a contract with the 
r Electromagnetic Radiation Laboratory, Air Force Cambridge Research Center. The Gaussian 
abscissas in Table II were computed by Murray Sherry using the ERD Applied Mathematics 
Group’s computer. 
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TABLE 1. Values of the Integrals 





























d feo das farota | seo abas | Sane 
0 | 3.6022 x 107 | 9.2361 x 10 | 3.0900 x 107 | 1.8134 

1 | 1.5318 X 10 | 2.9547 xX 10 | 1.7071 X 107 | 1.8013 

2 | 1.1372 X 107 | 1.6929 x 10 | 1.4551 x 10 | 2.9080 

3 | 1.1448 X 107 | 1.3333 x 107 | 1.9582 x 107 | 6.1454 

4 | 1.4091 X 107 | 1.3014 x 107 | 3.1323 x 107 | 1.5673 x 10 
5 | 2.0194 X 107 | 1.4915 X 107 | 5.7396 x 107 | 4.6231 x 10 
6 | 3.2683 X 107 | 1.9427 X 10 | 1.1815 1.5361 X 10° 
7 | 5.8546 X 107 | 2.8137 x 10? | 2.6921 5.6466 X 10° 
8 | 1.4433 4.4623 X 10° | 6.7053 2.2662 X 10° 
9 | 2.4148 7.6606 X 10? | 1.8073 X 10 | 9.8304 x 10° 
10 | 5.4566 1.4111 X 107 | 5.2283 X 10 | 4.5719 X 10 
11 | 1.3117 10 | 2.7692 x 107 | 1.6124 x 10° | 2.2648 x 10° 
12 | 3.3366 X10 | 5.7576 X 107 | 5.2719 x 10° | 4.8862 x 10° 
13 | 8.9412 10 | 1.2622 1.8189 X 10° | 6.5785 X 10° 
14 | 2.5146 X 10° | 2.9062 6.5962 X 10° | 3.8247 X 10° 
15 | 7.3980 x 10? | 7.0037 2.5077 X 10' | 2.3280 x 10° 
16 | 2.2704 X 10° | 1.7614 10 | 9.9395 X 10! | 1.4791 x 10° 
17 | 7.2503 X 10° | 4.6113 X10 | 4.1080 X 10° | 9.7850 X 10° 
18 | 2.4039 X 10' | 1.2537 x 10? | 1.7634 X 10° | 6.7240 x 10” 
19 | 8.2588 X 10‘ | 3.5333 X 10? | 7.8520 X 10° | 4.7898 x 10" 
20 | 2.9349 x 10° | 1.0302 X 10° | 3.3990 xX 10’ | 3.5306 X 10" 
ny fe? fa oo Sas | oo AF as fee Ra 
0 | 6.6995 x 10-1 | 1.5095 x 107 | 4.0149 x 107 | 2.1981 

1 | 2.9569 x 10-1 | 5.8235 x 10-? | 2.3627 x 107 | 1.8358 

2 | 2.0787 X 10 | 3.5690 X 10 | 2.2253 x 107 | 2.5158 

3 | 1.9019 X 107 | 2.8304 x 10° | 2.7450 X 10> | 4.6331 

4 | 2.1044 x 10 | 2.6958 X 10 | 4.1160 X 107 | 1.0577 x 10 
5 | 2.7098 X 10 | 2.9673 x 107 | 7.2041 X 10> | 2.8549 x 10 
6 | 3.9631 X 107 | 3.6859 x 107 | 1.4332 8.8270 X 10 
7 | 6.4691 X 107 | 5.0821 x 10% | 3.1793 3.0570 X 10? 
8 | 1.1649 7.6820 X 10-2 | 7.7504 1.6652 X 10° 
9 | 2.2771 1.2606 X 107 | 2.0531 4.8440 X 10° 
10 | 4.8106 2.2276 X 1071 | 5.8558 x 10 | 2.1679 x 10 
11 | 1.001110 | 4.2105 x 107 | 1.7851 X 10 | 1.0378 X 10° 
12 | 2.6334 10 | 8.4640 X 107 | 5.7806 X 10° | 5.2801 Xx 10° 
13 | 6.7332 X 10 | 1.8006 1.9785 X 10° | 2.8408 Xx 10° 
14 | 1.8155 X 10? | 4.0365 7.1263 X 10° | 1.6092 X 10? 
15 | 5.1394 X 10° | 9.4996 2.6914 X 10° | 9.5617 X 10? 
16 | 1.5244 X 10° | 2.3304 10 | 1.0624 X 10° | 5.9404 X 108 
17 | 4.7169 X 10° | 6.0109 X10 | 4.3707 X 10° | 3.8494 x 10° 
18 | 1.5193 X 10 | 1.6073 x 10 | 1.8714 X 10® | 2.5925 x 10” 
19 | 5.0821 X 10! | 4.4628 x 10? | 8.2973 X 10 | 1.8127 X 104 
20 | 1.7618 X 10° | 1.2840 x 10° | 3.8129 x 10" | 1.3127 X 10" 
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Transforming the first three integrals in (3) to the interval [—1, 1] we obtain: 


(4) 10) = 1 + Be + + [ 2Mz) ae, 
24 
in which 


{ 1 
If =o [ (1 + y)¥(40 + yl) dy, 


(5) ji" =4 [ (8 + 4{1 + yl)'¥(8 + 4[1 + yl) dy, 





N= { (16 + 4{1 + y)) (16 + 4{1 + y)) ay. 


\ 
The Gaussian quadrature formula with n abscissas, and weight function 1 gives 


(ef (1)): 
If = 4" 2 H,(1 + a;)*f(4{1 + aj), 


(6) I,” = 4 > H,<8 + 4{1 + a))*f(8 + 4{1 + a), 


Tis = 4 2) Hj(16 + All + aj))'f(16 + 4[1 + aj). 
j= 
Here the a; are the roots of the Legendre polynomial, P,(x), of order n, and 
the H; are the corresponding weights. These weights are given by: 


2 dP, 


“are Pe) = a. 


The integration procedure then is to select a value of n, obtain the values of the 
abscissas a; and weights H; , and evaluate the sums in (6). 





(7) H; 


3. Selection of n; Computation of the Integrals. Trial computations, using the 
tabulated values of the a; and H;[2] showed that a choice of n = 17 gives the 
integrals, for \ = 0, to five significant figures. For the higher values of X, it is neces- 
sary to take a larger number of intermediate abscissas in order to retain the five 
significant figure accuracy. On the basis of the following, admittedly naive con- 
siderations, we chose n = 27. 

First of all, the degree of precision of the Gaussian quadrature formula with n 
abscissas is 2n — 1, i.e. the approximation gives exact results if the integrand is a 
polynomial of degree 2n — 1 or less (On this point ef [1].). Since n = 17 gives 
“exact’’ five-figure accuracy for \ = 0 then f(4[1 + y]), in the interval [—1, 1] 
can be approximated to within a certain error by a polynomial of degree 33. Then 
the functions (1 + y)*-f(4[1 + y]) can be approximated by polynomials of degree 
33 + X. If these functions were actually polynomials of degree 33 + A, the number 
of abscissas necessary for exact integration would be given by 2n — 1 = 33 + A, 
or n = 17 + 4). Since, for our cases, A ranges from 0 to 20, we take n = 27. We 
have not been able to obtain an estimate of the error in using 27 points for \ > 0, 
but we feel from our experience with the computations, that the use of 27 abscissas 
retains the five figure accuracy of the results. 
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TABLE 2. Gaussian Weights and Abscissas 
n = 27. Weight function W(x) = 1. 








+a; Hj 

0 | . 114220849 
. 113972586 | . 113476364 
. 226459365 | . 111252467 
. 335993904 | . 107578263 
441148252 . 102501615 
.540551565 | .096088705 
632907972 | .088423136 
.717013550 .079604846 
.791771639 | .069748802 
. 856207908 | .058983516 
. 909482321 | .047449392 
. 950900558 .035297282 
.979923476 | .022686214 
. 996179263 | .009798976 





In order to obtain five significant figures in the integrals, it is necessary that the 
I,°, Is'*, Ti be correct to at least six significant figures. For this, because of round- 
off errors, it is again necessary that the H; and the values of f be correct to seven 
or eight significant figures, and this in turn means that the a; must be correct to 
eight or nine figures. Existing tables of the abscissas and weights do not contain 
the values with this accuracy. It is necessary then to compute these for our integra- 
tions. The a; are the roots of the Legendre polynomial 


13 . 
Px(z) = 2 (—1)"2™ q a ) ” 7 8) at 
and may be computed directly. The weights can be computed from (7). 

Table 2 shows the weights and abscissas for n = 27. The values are correct to 
nine decimal places. 

The computation of the function values f(4[1 + a,]), ete. were carried out 
using the British Association tables of the Airy functions [3]. In order to obtain 
these to eight figures it was necessary to use the second-order interpolation scheme 
suggested in the tables. 

As noted below the value of the last integral, 


[ 242) ae, 


in (4) is negligible, in our cases, in comparison with the sum J,° + Js + Ii. Con- 
sequently, the approximate integral is given by adding together the sums in the 
right side of (6). The results are shown in Table 1. 
4. Evaluation of [ a\f(x) dz. The contribution to the integral from the 
24 
‘tail’, 


[ a f(x) dx, 
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was estimated by integrating an asymptotic approximation to the integrand. This 
asymptotic estimate is derived from the asymptotic forms of the Airy functions 
given in [3]. In view of the fact that even the values of [7§ were negligible in all 
cases when compared with the sums J,” + J;'°, a satisfactory estimate of the integral 


: a f(x) dx 


is obtained by integration of the asymptotic form for f. The resulting integrals of 
these asymptotic expressions turn out to be incomplete gamma functions and 
may be evaluated with Pearson’s tables [4]. The result is that the ‘tails’ made no 


contribution to the fifth figure of the sum J,* + Io'* + %§, and are therefore neg- 
lected. 


Parke Mathematical Laboratories, Inc. 
Carlisle, Massachusetts 
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Recurrence Techniques for the Calculation 
of Bessel Functions 


By M. Goldstein and R. M. Thaler 


1. Introduction. The Bessel functions lend themselves most readily to calcu- 
lation by recurrence techniques [1]. Let us consider the regular and irregular Bessel 
function of real order and argument J,(z) and Y,(a). These functions both obey 
the same recurrence relation, viz. 


(1) P.a(z) + Fra(e) = = F,(2), 


where F(x) may be either J,(x) or Y,(x). If one is given Y,(x) and Y,4;(z) then 
Eq. (1) may be used to generate the functions Y,,,(2). For u > (x/2) the function 
Y,(x) increases extremely rapidly with increasing order, i.e., Y,(z) ~ (2u/zx)* 
and the functions Y,,,(z) calculated from Eq. (1) yield good accuracy for large n. 

However, if one is given J,(x) and J,4;(z), Eq. (1) gives poor accuracy for 
Jy4n(x), since for u > (2/2), Ju(x) ~ (2Qu/x)™. On the other hand, if one is given 
Jy4n(x) and Jyin4i(x), where n > (x/2), then one may again recur without loss of 
accuracy but this time in the direction of decreasing order. We shall first treat the 
problem of using the recurrence technique in the calculation of the regular Bessel 
function J,4,(z). Thus, let us find J,,,(z), forO S v < landn SN. 


2. The Regular Bessel Function. Consider a function F,,.(z), which obeys 
Eq. (1), and defined such that R 


Fy4uu(z) = 0 


P4u(Zt) = @ 


(2) 


where a may be chosen to be any constant, and M > N. By successive application 
of the recurrence relation, Eq. (1), we may now generate F,,4i(z), --: , F,(2). 

Since F,4.44:(2) and F,,4(z) can be treated as the same linear combination of 
the regular and irregular Bessel functions, then 


Py ugi(t) = advamgi(2) + BY r4.041(2), 
Fy4u(x) = od 4. u( 2) + BY,4u(2), 


(3) 


and, in general, 


(4) Fy 4n(2) = ad y4n() - BY y4n(X); 
so that 
(5) Py 4u(x) = od »4.0(2) [1 + g Foawle) |, 





Received Jan. 12, 1959. Work performed under the auspices of the U. S. Atomic Energy 
Commission. 
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CALCULATION OF BESSEL FUNCTIONS 














and 
F,4x(z) ™ oe] y4.0 (2) [1 + B pete) | 
a J 4"(z) 
Since M is chosen such that M > N, then it is clear that 
B Y,46(z) 
J +4.0(2) Yogw(x)\ (J r4a(x) 
(6) AL. ( + <1, 
B Yi4a(z) Y,4m(z) J +4.(Z) 
a Jr4u(z) 
so that forn <= N 
(7) Fyan(2) SS ad r4n(2). 


Clearly M may be chosen so large that F,,,(z) = aJ,4.(z) to any desired nu- 
merical accuracy. 

Determination of a will then yield the regular Bessel function J,,,(z) for 
n & N. To do this one may make use of one of several addition theorems; for ex- 
ample, the addition theorem [2]: 


(8) 2’ > (v + 2m) J42m(xz)x°T(v + m)/m! = 1. 
For v = 0 Eq. (8) reduces to the familiar result 
(9) Jo(xz) +2 2) Jom(x) = 1. 


To any given accuracy there exists an even integer L, such that 


L/2 


(0) 2” > (vp + 2m) Jrsam(z)2-°T(» + m)/m! & 1. 
m=0 


If L S N, then to the desired accuracy one may write 
L/2 


(11) 2’ > (v + 2m) Fys2m(x)z°T(v + m)/m! & a. 


m=0 


If L > N, the desired accuracy may nevertheless be obtained by increasing M. 
Eq. (11) may be rewritten in a form more suitable for numerical computation as: 


L/2 


(12) 2, mF 42m(2) = a, 
where 


7~— (2) ra #05 
Zz 

_ (vy + 2m)(v + m — 1) 

m(v + 2m — 2) 


Om 





Pm-1 . 
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3. The Irregular Bessel Function. The irregular Bessel function Y,(zx) is defined 
as 


J,(x) cos vr — J—(t) 


sin vr 


(13) Y,(z) = 





Expansions for J_,(x) in terms of J,,m(x) are readily obtained [2], for example 
one may write: 


J_.(z) = (5) T(1 — 2v) x (y ot) 
(14) 
I'(v + m) 1 


‘Ta —m — 2v) T(1 + m — ») 
Substitution of Eq. (14) into Eq. (13) yields the result that 





J y42m(2). 


(15) Y,(z) = 2. Ym J v42m( 2), 
where 


1 2) r’(1 + ») 
cot yr — —|-} ————, 
TT 





-m 
zt v 
m= ()(@) (FE) rato. 
mw} \x 1 
“ _(v + 2m)(2v + m — 1)(v + m — 1) 
ms m(m — v)(v + 2m — 2) _— 


For smal] values of | v | the coefficient ye may be expandéd as 


‘Se . 2 a 2A° 
Yo = =| 4 — “(+ A’ +t )+y (3 + AS, + *) 





(16) a 
= f (H+ Sy AS, + vy 7) + |, 
where 
A = 0.577 215 6649 --- + log 5 
and 


p=1 
x 

So = S’ 

S; = 1 .202 056 903 ---, 
= 

S= 5: : 
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Thus, if we have obtained J,,,(z), we may use Eq. (15) to caleulate Y,(z). 
The Wronskian relation, 


i) 


(17) Y,(2)J 4a(z) — Youa(z)J.( 2) = —, 


us 


& | 


gives Y,4:(2). From Y,(x), Y,4:(2) one may obtain the values of Y,,,(z), n > 1, 
by recurrence, Eq. (1). 

If x is close to a zero of J,(xz), then Eq. (17) does not provide a suitable method 
for obtaining Y,4;(2). To avoid this difficulty one may obtain a series for Y,4,(27) 
which is analogous to the series for Y,(x), Eq. (15). This is readily accomplished 
by differentiating Eq. (15) and using the relation: 


(18) Visa) =” ¥,(x) — 1 (2), 
x dx 
After some manipulation one may obtain the result 
(19) Yess = 2D bdese, 
m=0 
where 
1 2 1+2» 
w= -1(2) r(1+y), 
w\z 
1 
i = vo — 5 7Yi> 
bom - = Ym > m > 1 ? 
2 


1 
Eamsa = 5 (Ym — Ymsi) » m2 I, 
and the ym are as defined by Eq. (15). 

4. Bessel Functions of Imaginary Argument. Analogous formulae are easily 
derived for the Bessel functions of imaginary argument J,(x) and K,(x). The regular 
function J,(2) obeys the recurrence relation 

2p 
(20) I,1(z) +a Tyai(z) = m9 I(x) . 
The irregular function K,(x) obeys the relation 


(21) K,a(2) — Kra(2) = — 2 K(2). 


It is convenient to define the function K,(x) = (—1)’K,(x). Then J,(x) and K,(x) 
both obey the same relation, viz: 


(22) G,a(2) — Gryalz) = ~ G(z) 
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One readily sees then that taking 
O = Grows = olyymai(t) + BRoyms(2) 
@ = Gryu = aly4u(z) + BR, 4u(2) 

will once again yield 


(24) Gy4n(2) = al,sn(Z), 


(23) 


to any desired accuracy for M > N 2 n. The addition theorem analogous to Eq. 
(8) is 


(25) 22 (-1)"( + 2m)I,42m(x)x°T(v + m)/m! = 1. 


In order to avoid the use of an alternating series, however, it proves useful to use a 
different addition theorem [2], viz: 


2VT(i+v) S (v+m) ~z . 
(26) 2 (2) Ti+ 2) yy = T'(2v + m)(e*D,4m(z)) = 1. 


This leads to the following formula for a: 





(27) a SQ VnGrim(2), 
where 
hme (2) r(1 + ») 
x 


D (vy + m)(2v + m — 1) 
m(v + m — 1) 


Vmn—1 . 





Yen 


The irregular function K,(x) may be treated completely analogously to Y,(z). 
The irregular function K,(z) is defined as 


(28) K(2) = [FH = Ae) |. 


2 sin yr 





By means of the expansion 





2 m! I(l— m — 2p) 
(29) , 
. Ta +m—») T,42m(x) ’ 
one obtains the result that 
(30) K,(x) = 2) dmlrs2m(2); 
where 


6b Oro] 
vy _SIN vr zx 








Eq. 


se a 


x), 
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(vy + 2) 
s=FUr oy i= 


(v + 2m) (2 + m — 1)(v + m — 1) 
m(v + 2m — 2)(m — v) 





bm = 





5m4t - 


For small values of | » | the coefficient 5) may be expanded as 


. a A* 
& = —|A—»(2 + A? -* +9 + AS + 72 


S 2A8S S 7x 
= (4% ese 
where A and S,, are defined as in Eq. (16). 

Thus, as before, if we have obtained /,,,(2) we may use Eq. (30) to calculate 
K,(z). 

The Wronskian relation, 


(31) 





(32) K,(x)I,4:(z) + Kys:(x)1,(2) = 


gives K,4:(2). From K,(x), K,4:(2) one may obtain the values of K,,(x) for 
n > 1 by recurrence Eq. (21). Unlike J, , J,(x) is not an oscillatory function and 
the Wronskian can, therefore, serve to give K,,;(z) without difficulty. 

However, Eq. (30) does not yield high accuracy for x > 1, since for large x 


(33) bol, ~ — ai Se 


This difficulty can, in practice, be overcome by the evaluation [3], [4] of the integral 
representation of K,(zx) 


(34) K,(z) = | e =" cosh vy dy. 
0 


5. Large Values of the Argument. If the values and derivatives of any one of 
the functions J,(z), Y,(x), I,(z), K,() are readily obtained, then the above 
techniques are somewhat too cumbersome. In particular, if one has Y,(z), 
Y,4:(x)[K,(z), K,4:(x)] then one may obtain the values of Y,+.(2)[K,4.(x)] by a 
straightforward recurrence. On the other hand, if one has J,(z)[J,(x)], then one 
may follow the procedure of Eqs. (2-7) [Eqs. (23-24)]. However, now a is given 
simply by 


(35) a = F,(x)/J,(z) 
for J,(x), or 

(36) a = G,(x)/I,(z) 
for I,(x). 


For x 2 10 the necessary values and derivatives of the functions J,(xz), Y,(z), 
I,(x), K,(x) are easily obtained by the phase-amplitude method [5] for 0 S v < 1. 
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For very large values of the argument, this technique is more suitable for computa- 
tions than the methods outlined in the previous sections. 

A subroutine [6] for a high speed calculating machine, the IBM 704, has been 
written incorporating the methods described here for the calculation of Bessel 
functions. 


New York University, New York City, New York and Los Alamos Scientific Laboratory, Los 
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Calculation of Transient Excitation of Ship Hulls 
by Finite Difference Methods 


By Harry Polachek 


1. Introduction. A system of finite difference equations based on the non-uniform 
beam theory is developed for use in the calculation of the response of a ship hull 
to transient forces. The conditions for stability of these equations (and hence the 
conditions for validity of the numerical results) are derived. The feasibility of the 
method is tested by the solution of a vibration problem for a specific hull, the 
details of which are discussed in [1] and [2]. The use of this method lends itself to 
the solution of a wide class of problems related to the structural design of vessel 
hulls or other structures subject to transient forces. The solution has been pro- 
grammed and carried out on the Bureau of Ships UNIVAC System, Applied 
Mathematics Laboratory, David Taylor Model Basin. 


2. Governing Equations. The equations governing the motion of a ship hull 
based on uniform and non-uniform beam theory as developed by Timoshenko [3] 
and others are discussed in considerable detail by R. T. McGoldrick and V. L. 
Russo in [4]. A finite difference method for obtaining numerical solutions to these 
equations is presented here. For this purpose, the equations describing the damped 
vertical (or torsion-free horizontal) excitation of a ship hull subjected to a transient 
force will be used. The results may be directly extended to more general types of 
motion. The system of partial differential equations describing this type of motion, 
as given in [4] is: 





ay  ! 
(1) Uae tea te = P(z,t) 
ay » Oe . 
(2) | ra +1 _ = 0 
oy 
(3) M = (EI) — 
Ox 
“ - oy 
(4) V = (KAG) y — (KAG) = 
Ox 
where, 
t = time 
x = distance coordinate along the longitudinal axis of a vessel 
y = displacement normal to the longitudinal axis 
Y = rotation of transverse section about an axis normal to the (zy) 
plane (z-axis) 
M = bending moment 
4 = net shear force in y direction 
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m = apparent mass (per unit length) 

c = damping factor (per unit length) 

P = force (per unit length) acting upon ship hull 

Ins = mass moment of inertia about z axis (per unit length)—‘“rotary 
inertia” 

(EI) = bending rigidity factor 

(KAG) = shearing rigidity factor. 


3. Finite Difference Representation. The set of equations (1) to (4) constitutesa 
system of partial differential equations which is assumed to govern the motion of 
the hull of a vessel, as simulated by a freely vibrating non-uniform beam. If the 
state of motion of the ship hull is known at any time t , it is possible by obtaining 
the solution to these equations to determine its motion at any subsequent time. 
This system of equations may be represented approximately in finite difference 
form by replacing the partial derivatives by equivalent ratios of small finite incre- 
ments. In making this substitution we will use the following notation: 

tc——initial time 
Xo initial position 
At = increment in time 
Ax = increment in length 
Ay, Ay, AM, AV = corresponding increments in the dependent variables y, y, 
M and VJ, respectively 
In = % + ndAz 
Yn+a = ylxo + (n + 4) Az, to + mAt] = value of y at the position 
ao + (n + 4) Az and at the time & + mars” ul : 3): 
M7" = M[xo + nAz, to + (m + 1)At] = value of M at the position 





. n=0,1,2-:-- 
xo + nAz and at the time t + (m+ 1)avs{ on 1.9--: 


Similar notation will be used for other variables and factors in the equations. For 
quantities which do not vary with time, the superscript will be omitted. 
y 





if eo aoe 





















































Fig. 1a.—Ship hull. 


Ko KL X) Mt %20 










































































Fie. 1b.—Representation by non-uniform segmented beam. 
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Figures (la) and (1b) show the representation of a ship hull by means of a 
non-uniform beam and its division in twenty increments, as used in the initial 
solution. Using the above notation, the displacement of the hull at the position z. 
(the boundary point between the second and third increments) at the end of the 
mth time increment, for instance, will be represented by y2”; its displacement at the 
center of the second increment will be represented by yi, . The average mass per 
unit length at the center of the twentieth increment will similarly be represented 
by p19.5 - 

To solve the above system of equations we propose to make the following specific 
finite difference substitutions in the system of equations (1) to (4): 











(5) (z2)’ _ ath — Ques + 9 
OP / n44 (Az)? 

(6) (2) _ vat} — Qyne + yea 
OF J nti (At)? 

(7) (2) _ nea — yea 
oV 4 — Ve+t a y. 

” (2) * Az 

(9) (ary’ _ Min-M." 
Ox J nti Ar 
dy\"" —oyntt — ytti 

- ( uy’ we Ar 

o (2) = 


The resulting system of finite difference equations is given below: 


™m m m— Cn m m— 
ynth = (2yn44 — nat) 7” (=) (yn+4 — yns4) At 
n+ 











(12) ™m m 2 2 
= (Viw— Van ) (At) + pe ( At)” 
Mn+} Ar ott Mn+} 
m+1 __ m — awl M41 — M,” (At)? — Vin+y( At)” 
(13) Yn+4 ag (Qyr44 Yn+4) + (Tye) n+ Az (Tye) nat 
m+1 m+1 
(14) Mz = (BI), 22 te 
Az 
m+1 - m+1 
(15) Vi"! = (KAG), 73" — (KAG), 244 —¥et 
where 


Yn” = (nua + Yn), Vasa = 3(V2" + View). 


In order to define completely the motion of the ship hull (non-uniform beam) it 
is required that initial and boundary conditions be fixed. If we begin our computa- 
tion when the vessel is at rest (just prior to subjecting it to any force) we have the 
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condition at & = 0, y = y = dy/dt = dy/dt = V = M = 0. At all time, ¢, we 
have the condition Mp = Vo = Ma = Va = 0, on the basis of our assumption 
that the hull is vibrating freely. 

The above system of equations (equations (12) to (15)) taken together with 
the supplementary initial and boundary conditions may be summarized in the 
following form used in programming the solution on a high-speed digital calculator 
(UNIVAC System) : 

(12’) yety = (1 + Ki)(ynay — yaad) + you + Ka(Viu — Vin") + Ks 

(13’) yede = (nay — Yash) + ea + Ks(Mia — Mn”) + Ke( Vi + Vn") 
(14) MR" = Ka(yndy — yn44) 

(15°) Va*t = Kr(yety + yet) + Ke(ynti — ys44) 

where K, , K2, --- Ks are multiplying factors which are, with the exception of K;, 
functions only of the characteristics of the vessel hull and the selected time and 
space intervals, and which may be precomputed prior to the main calculation. 
K; may also be precomputed in the case of the application of a constant force; in 


the case the acting force is a function of time a new set of values K; must be com- 
puted for each time interval. Specifically, 


la 


(At)? Pon 














a 2 a ee » K= (af), 
Mn+} Mn+} Ax Mn+ 
(EI), 1 (At)? (At) 
16 K, = ——_, , a ee K.= — : 
( ) | aes Ax . (Iyz) n+ Ax P 2(L uz) n+ 
_ (KAG), _ _ (KAG)s ; 
at. *. =e Toe - 





\ 


The initial and boundary conditions are given by the relations, 


Yaa = Yet = ¥044 = eH = Va = M, = 0, 
(17) 
| Mo” = Moo = = 730 = 0. 


From the above it is also possible to derive the following useful relations for the 
end-point values of y: 


(18) yo" = ”" 20 = Yi9.53 Yor = ys = 3 "5 y20 = Yio + S vis. 

4. Numerical Stability. In order to carry out successfully the solution of a system 
of partial differential equations such as equations (1) to (4) by finite difference 
methods, the finite difference equivalent system (12) to (15) must be stable in the 
sense discussed in [6], [7] and [8]. We will now derive the conditions under which this 
system of equations will satisfy these stability requirements. The conditions of sta- 
bility are satisfied if the amplitude of a small disturbance, introduced at any time, ¢, 
does not increase exponentially with successive time steps. This condition may be 
stated as follows: If 6F (2, t) and 6F(a, t + At) are values of a variation (or per- 
turbation) of any of the dependent variables y, y, M and V in the system, then it is 
said to be stable provided | F(a, t + At)/éF(z, t) | S |. To determine the condi- 
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tions for stability we will introduce perturbations dy, dy, 5M and éV in the de- 
pendent variables y, y, M and V, respectively. Substituting in equations (12) to 
(15) we obtain the variational equation system: 


pins4(5ynt4 — 26y44 + Synz4) + Cnsy(Syn4y — Synz4) Af 
(19) 








+ a (8V241 — 8V.") = 0 
I n m ™m m— 6M; = 6M,” 7m 
(20) “as (Synth seald 267 44 + by'n34) _ or “hn + bV aH = 0 
(21) at (dynti — byn44) — 6Mz™ = 0 
by nid oe Syn} oe m+1 sve" i 
(22) aa yn + (KAG), e 


We will assume in this analysis that within a small region in the (z, ¢) plane the 
coefficients (yz, I,., etc.) of the variational functions can be treated as constants. 
A solution of the system of equations (19) to (22) can then be obtained in the form 


yn” = ae Pt temss 
wa » ~ be'Patamse 
(23) ie Pit eg Pa tamse 


m iBn+amAt 
— «<6 


where a, b, c, d are real constants and a complex. Substituting (23) in (19) to (22) 
we obtain a system of linear homogeneous equations for the quantities a, b, c and 
d which has a non-trivial solution provided the determinant D of the coefficients is 
identically zero, where 








| w(A— 242-2) + cat(1 — 2-1) 0 0 um les —e-4) 
Az 
, a 
0 —.. (A—2+ 7) _ (3 — ce?) 1 
(at)? Az 
(24) D = oe 
0 — (e2 — 2) =] 0 
Az 
ig 8 
(er —e-F a 0 Pé 
Az (KAG) 
and where \ = e““‘. From the above we obtain 
- “as EI ( ais a) 
~1 
(25) _ Te(A — 2+ ] 
(KAG) (At)? 





— (20% (ain 8 [ EE (a in 8) — 12 = 249) 
(Az)! (2i sin 5) | ( £L (2 sin *) (an? 0. 
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If we further assume At < (u/c) we obtain, 


(25’) (eh, ue + 4r sin +B) ar(E1) sin’ — 8 5 tLe :) + pt(At)? = 0, 


KAG 
or 
28 u(ET) 2 2 . B_ 
(26) pe + [a r sin’ 8 (tn + ED) + w(as? |] s+ 16°CED) sint 8 = 0 
where 
a Gh -1 _ (at)? 
—=AX 2+,” , = (a) 


Equation (26) above is a fourth degree polynomial equation in \. If we designate 
one of its roots \; we draw the conclusion from the relation § = \ — 2 +X” 
that A» = 1/A;, will be another root. For stability both |A,| S land || = 


|1/.| S$ 1. It follows that || = |A2| = 1. Likewise, | A3| = | | = 
Let A. = cosy +7 sin y; Ax = cosy — i siny = 1/d,, then, 
(27) & = 2(cosy — 1), 
or 
—-4Sé<0. 


On the other hand, from equation (26) we obtain 


(28) = 548 SUS + w(at)*) + (UT + 2u(a1)*US + 2(At)*)4 


and 
(29) ¢= si [(US + w(At)*?) — (U°T? + Qu(at)?US + y?(at)*)4) 
where, 








It follows that, for stability, the inequalities 











(30) -4 5 -Ea6 + 2u(At)?US + w2(At)*)] < 0 
and 
(31) -4< 3 (US + w(at)*) — (U*T* + 2u(at)?US + 2(at)")) < 0 


must hold. It may be seen by examining the expression for £ that its value is always 
less than or at most equal to zero and that the left hand inequalities (30) and (31) 
are satisfied if 


(32) (US + u(at)*) $4 re <2. 


1 








+ 


WA 
o 


<0 


lways 
(31) 
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In addition, from the previous discussion it follows that At must be chosen so that 
(33) At « = , 


The above result may be stated as follows: 


THEOREM 


The finite difference system, consisting of equations (12) to (15), governing the 
transient motion of a ship hull, is numerically stable provided the time increment 
At is chosen sufficiently small so that At < w/c and 


(32) (at)* Hl (A) 
~ Ty(KAG) + w(EI) + .25u(Ar)*(KAG) 
throughout the range of solution. 

Numerical instability is usually accompanied by violent variations in the com- 
puted functions which invalidate the solution. This effect is shown in Figure 2. 
Whereas a time interval satisfying inequalities (32) and (33) guarantees stability, 
it should be pointed out that the interval need not necessarily be as small as indi- 
cated by these inequalities. 





5. Computational Procedure. We are now prepared to use the system of finite 
difference equations (12’) to (15’) and (16) to (18) to calculate the motion of a 
ship hull in response to a force acting upon it. The physical characteristics of the 
ship hull and the magnitude of the forces acting upon it (ie. u, C, P, Iy., (EI) 
and (KAG)) are calculated on the basis of the theory of elasticity from experi- 
mentally determined physical quantities. (For instance, see Table 7, [5] for values 








0.002 
Violent Variation due to 
Numerical Instability 
(AT =0.012) ey 
0.00! N 
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Fig. 2.—Numerical instability resulting from incorrect choice of integration interval. 
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of these parameters in the case of the SS Gopher Mariner.) From these the maxi- 
mum value that can be assigned to At, the time increment, in order to insure nu- 
merical stability may be calculated on the basis of the inequalities (32’) and (33). 
Then the coefficients K,, K.,--- Ks are calculated by use of equations (16). 
These should be listed at the full-interval positions or at the half-interval positions 
(See Fig. 1b) as follows: 


K, , K2, Ks, Ks , Ke——iisted at half-interval positions 
Ky, Kx, Ks listed at full (integral) interval positions. 





From the known conditions (y, y, M, and V) at time ¢ = 0 (equation 17) we 
proceed to calculate the values of y and y at t = At, by use of equations (12’) and 
(13’). We then obtain the values of M and V at t = At from equations (14’) and 
(15’), also using the boundary relations My = Mx = Vo = Va = 0. We may 
now repeat this cycle any number of times, obtaining the values of the variables 
y, y, M and V att = 2At, 3At, --- , etc.—until we reach any desired value of time, t. 

The above computation procedure was programmed for solution on the UNIVAC 
system, and trial solutions carried out. The results appear successful in every re- 
spect. Calculations at varying time intervals demonstrate the feasibility of pro- 
ducing accurate solutions much beyond engineering requirements. 


Applied Mathematics Laboratory 
David Taylor Model Basin, 
Washington, District of Columbia 
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TECHNICAL NOTES AND SHORT PAPERS 


Skip-Term Summation of Sequences 
By Irwin Roman 


The literature on the summation of sequences and the finding of the sums of the 
corresponding series is voluminous and well known. For rapidly convergent se- 
quences term-by-term evaluations furnish no great obstacle. Likewise, if the se- 
quence can be matched with a sequence whose sum is known and the term-by- 
term difference of the two sequences converges rapidly, this difference-sequence 
can be evaluated by terms. 

For series that converge slowly, various devices have been suggested. If the in- 
tegral and derivatives of the general term, expressed as a function of the term index, 
can be evaluated, the Euler [3]-Maclaurin [7] formula or the extended-range form 
of these formulas [9] can be used to find the sum of the series, or of the residue 
after a selected term. 

If a series converges sufficiently smoothly, the formulas of Lubbock [6], as mod- 
ified by de Morgan [2], can be used to evaluate a partial sum of the terms of the 
series from those terms which are separated by a selected gap. Lubbock’s formulas 
were revived by Sprague [12], have been discussed in Whittaker and Robinson [14], 
and were extended by Steffensen [13]. Lubbock’s formulas all involve differences, 
and tables have been published by Davis [1] and others, to facilitate the evalua- 
tions. 

The limit of the sum of a slowly convergent infinite series has been approxi- 
mated by Salzer [10], who selected the reciprocal of the sequence index as the inde- 
pendent variable in Lagrangian interpolation and found the approximate limit 
as this variable approaches zero. He selects the last m members of the sequence of 
partial sums from the n members which have been calculated. Simple factors were 
published for n = 5, 10, 15, and 20, and for m = 4, 7, and 11, where m < n. 
Later, Salzer [11] published formulas and coefficients for the sum of a finite series, 
and Horgan [4] prepared decimal tables of coefficients for finite sums. These for- 
mulas involve no differences and are convenient for many purposes. However, they 
involve extrapolation. The present formulas involve only interpolation. Only a few 
coefficients are included, as examples of the method. 

If a sequence of terms {} can be approximated by a polynomial of degree n 
in the index k, the values of (n + 1) arbitrarily selected terms of the sequence 
can be used to determine the value of every other term. If {y,} is a subset of {}, 
containing (n +1) elements, Lagrangian interpolation determines 


ae Dy Au"(k) Ys ’ 


where A,,"(k) can be determined independently of the values of y,. The sum of 
the sequence of m terms of {m} is S = >-him = Doha Dot-o Au"(k)yu. As 
the number of terms is finite, the order of summation can be reversed, so that 
S = >}°2_. B."y., where B,” = >-%, A."(k), can be determined independently 
of y, . 





Received June 28, 1957; in revised form March 18, 1958. 
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Although the preceding analysis is valid for every subset of the set, an applica- B 


tion of special interest is that of skip-term summation in which the subset {y,} 
is determined by the relation y,, = ,.. The first term of the set and of the subset 
is Yo = m and the final term is y, = m2. For this choice m = gn and the terms 
of the subset are taken from those of the set with a gap g in the index. The initial 
term is not included in the sum, but its value is needed to obtain the sum. This 
choice permits the summation of a sequence of terms in blocks, without including 
overlapping terms twice, and is convenient in evaluating the sum of a smooth 
sequence after selected terms have been computed and added. For this choice, 
the sum of the ng terms of the original sequence is S = ) 2» Bu". , where 
B,” = 2 3°, A.(k). 

The coefficient A,,”"(k) is a polynomial of degree n in k, and depends on u. Accord- 
ingly, the summation involves sums of the type ) i., k’, which can be shown 
[9] to have the value 


t w be 1 v—2u+l 1 v—2u+l 
_ u . - we i. 
2B = ee =o FI! l(: - 5) (5) |. 


where w is the largest integer in v/2 and where by = 1, bs = —4, bs = g%q, and 





bow = — >, bou-o:/(2i +1)! foru > 0. 
i=l 
The sum can be written 


t v+l w 
Lk = (5) Do deulv!/(v — 2u + 1) !Nelas2u, 
k=1 u=0 


where g;' = (2+ 1)’ — 1. Specifically, for the sum from 
k=1 to k=t, Dk=@/8 DOF = (e: — o)/2%, 
Dok = (4 — 22)/64, DE k* = (395 — 10gs + 7¢1)/480, 
Lo = (es — 5, + 72) /384, 
Li k° = (397 — 2les + 4993 — 31¢:)/2688, 
Dd k" = (39s — 28y5 + 98e, — 124g) /6144, 
Xk = (Seo — 60g; + 29495 — 620y; + 381y1)/23040, 
dX & = (¢1 — 15¢s + 98¢s — 310~, +381¢2)/10240, 
> hk = (391 — 55¢9 + 4629; — 20469; + 41919; —2555¢; ) /67584, 
where the value of ¢ is implied. The customary forms of these sums are polynomials 
or products of polynomials in t. (See, e.g., Jolley [5] formulas 15 to 22 inclusive.) 
The Lagrangian coefficients can be obtained in various manners. One method is 
that of expanding into polynomials the customary forms involving products. 
Another method involves the assumption of a polynomial of the proper degree and 
the determination of the coefficients by solving the system of linear equations cor- 


responding to the various values of the terms of the subset used in the summation. 
Forn = 2p, the M.T.P. tables [8] of Lagrangian interpolation coefficients are given 
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for the range u = [—p(1)p] so that the terms of the original set may be taken on 
the range k = [—gp(1)gp]). For this choice k = gu. The subset is {y,} and, drop- 


ping the superscript n, the desired sum is 


oP P oP 
S= 2 > x Buyu where B, = p> A,(k/g). 
=l-op 


For n even, A_.(—z) = A,(z) so that 
B, = Au(0) + A(p) + 2 [A-a(k/9) + Au(k/9)] 


For integral values of uandv, A.(u) = land A,(v) = Oforv + u. Hence 


op—l 


B_, = 2 [A-a(k/g) + Au(k/g)], 
B,=1+ B_,, 
By = 1+ 22 Ad(k/9), 
op—i 
B_, = B, = 2 [A-u(k/9) + Au(k/g)]) for 1Susp-—1. 


Accordingly, the sum of the original set is 
p-l 


S = i + Boyo + By + 2d, Bu(ya T Yu). 


VaLuEs or B, 
Three Point 












































\ ‘ 5 | 10 | 25 | 50 100 
a 
oa 1.2 | 2.85 7.84 16.17 | 32.835 
0 6.6 | 13.30 33.32 66.66 | 133.330 
+1 2.2 3.85 8.84 17.17 33.835 
Seven Point 
“A 5 10 25 50 100 
—3 {1.007 808 | 2.450 365 875 | 6.830 150 082 56 |14.147 218 165 O08 | 28.787 894 830 158 75 
+2 |7.594 752 |15.368 629 750 |38.547 432 304 64 |77.130 857 609 52 |154.279 714 344 047 50 
+1 |1.150 720 | 2.022 188 125 | 4.858 920 038 40 | 9.661 606 076 20 | 19.295 O89 152 381 25 
0 |9.493 440 |19.317 632 500 |48.526 995 148 80 |97.120 636 298 40 |194.274 603 346 825 00 
+3 |2.007 808 | 3.450 365 875 | 7.830 159 082 56 |15.147 218 165 08 | 29.787 894 830 158 75 
Eleven Point 
\——__}—_. 5 10 
0 +34.028 105 216 +70.538 954 406 359 375 00 
—5 +0.892 014 592 +2.208 634 246 136 718 75 
+5 +1.892 014 592 +3.208 634 246 136 718 75 
+1 — 20.316 441 600 —42.788 478 982 539 062 50 
+2 +21.805 521 920 +45.019 765 486 718 750 00 
+3 —3.588 049 920 —7.870 561 162 597 656 25 
+4 +8.692 902 400 +17.661 163 209 101 562 50 
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Values of B, are shown for several values of n and g in the accompanying Table. 
The values are exact and have been checked by the relation > >?__, B, = 29p 
based on y, = 1. With allowance for rounding errors in the tenth decimal place, 
they have also been verified by addition of A; from M.T.P. 

If m ¥ gn, the method can be used only when the sequence can be divided into 
blocks each of the selected form. The sum of the terms ahead of the first block must 
be calculated separately. For example, for sixth order summation of the block 
k = [—39(1)3g], the sum from 1 — 3g to 3g is 


S = Bsy-s + Boyo + Boys + Bi(ya + m1) + Bo(y-e + y). 
For two consecutive blocks k = [0(1)12g], the sum from 1 to 12g is 
S = Bs(yo + ys) + Bolys + yo) + Bays + yr2) 
+ Bily2 + yo t+ ys + yro) + Boys + ys + yx + yu). 


For n = 2p — 1, the M.T-.P. tables [8] use the subset range u = [1 — p(1)pl, 
so that the original set must be taken on the range 


k = [g{1 — p} + 1(1)gp). 
For this case, A;_.(1 — x) = A,(x). Then 
gp op 
S= > -> Buy. where Bu= )>. A,(k/g). 
k=1— tie u=0 k=1—(p—1)9 
The value of B,, can be written 
gp oPp—9 Bhilai 
Bu= DO Autk/g) + D [Aulk/g) + Aug — 1+ &/g)).- 
k=gp—g+l k=1 


Accordingly the method is somewhat more involved for n odd than for n even. 

Tests of this method on some of the examples given by Salzer [11] indicate a 
superiority of the latter method in summing slowly convergent series. However, 
the present method is convenient in summing a finite sequence of terms that are 
sufficiently smooth but the infinite series of which does not converge. It is also 
convenient when individual terms are not computed readily except for multiples 
of a constant, and where sequences of terms are simpler to compute for large index 
than for small. 


United States Geological Survey 
Washington, District of Columbia 
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On a Conjecture Concerning the Primes 


By R. B. Killgrove and K. E. Ralston 


Consider the sequence { Po;}, i= = 0, 1, 2,--- 
ber, Poo = 2, Pa = 3, Pe = 5 ov; 
primes by the recursion relation 

Pi; = | Pin — Pins]. 
The conjecture (Norman L. Gilbreath, private communication, July 1958) is 


then that Pi = 1 for all i > 0. The validity of the conjecture for the first few 
primes can be seen from the following table of their absolute differences. 


p> Se PORE ules 


, where Po; is the jth prime num- 
. Now duties the absolute differences of the 


a 2 2 
ape © 0 
i @ 0 
;' 2 
1 


There are an uncountable number of sequences {bo;} with the property that 
their absolute differences bi defined as above are unity. In particular the sequences 
{k + 1, k, k, --- } and any sequence of the form {bo = 1; bo; = 0 or 2,7 > 0} 
have this property. Furthermore it can easily be verified that any sequence, {bp;}, 
with the required property has its first absolute differences bounded by the se- 
quence {2%}, that is, bi; < 2’. 

Consider again the absolute differences of the primes. Since all primes greater 
than 2 are odd numbers it follows that all differences P;; , 7 > 0, are even numbers. 
Now, if for some 7 and all 7, 0 < 7 < M, we have P;; = 0 or 2 and Pj» = 1, then 
all of the differences that derive from them will be bounded by 2, from which it 
follows that 

Pio, Pisio, Piseo, +++» Pitmayo = 1. 


We now define the function P(i) to be the largest integer M such that P;; < 2 
for all 7 < M. Thus we can say that Pi = 1 fori S k < P(i) +7. 

A routine was coded for the SWAC to evaluate this function P(7), using the 
primes less than 792,722 from a sieve prepared by D. H. Lehmer. The results of 





Received Oct. 7, 1958. The preparation of this paper was sponsored by the Office of Naval 
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this calculation are shown in the following table. From these results it is seen that 
the conjecture holds for all primes less than 792,722, which amounts to the first 











63,419 primes. 
TABLE OF THE FuncTION P(i) For 0 < i < 95 
i PG) | PG) +i i Pi) PGi) +i i PG) | Pii) +i 
1 3 4 33 867 900 65 23266 | 23331 
2 8 10 | 34 866 900 | 66 23265 | 23331 
3 14 17 35 2180 2215 67 23264 | 23331 
4 14 18 36 2179 2215 68 23263 | 23331 
5 25 30 37 2178 2215 69 31500 | 31569 
6 24 30 38 2177 2215 70 | 31499 | 31569 
7 23 30 39 2771 2810 71 31498 | 31569 
8 22 30 40 2770 2810 72 31497 | 31569 
9 25 34 41 2769 | 2810 73 31528 | 31601 
10 59 69 42 2768 | 2810 74 31527 31601 
11 98 109 43 2767 2810 75 31526 | 31601 
12 97 109 44 2766 2810 76 31526 | 31602 
13 98 111 45 2765 2810 77 31528 | 31605 
14 97 111 46 2764 2810 78 31527 | 31605 
15 174 189 47 2763 2810 79 31526 | 31605 
16 176 192 48 2763 2811 80 31526 31606 
17 176 193 49 2763 2812 81 31536 31617 
18 176 194 50 2763 2813 82 31535 31617 
19 176 195 51 3366 3417 83 31534 | 31617 
20 291 311 52 4208 4260 84 31533 | 31617 
21 290 311 53 4207 4260 85 31532 | 31617 
22 289 311 54 4206 4260 86 31531 | 31617 
23 7 763 55 4205 4260 87 31538 | 31625 
24 874 898 56 4204 4260 88 31537 | 31625 
25 873 898 57 5943 6000 89 31536 | 31625 
26 872 898 58 5944 6002 90 31535 | 31625 
27 873 900 59 5943 6002 91 31534 | 31625 
28 872 900 60 5942 6002 92 31535 31627 
29 871 900 61 5941 6002 93 31534 31627 
30 870 900 62 5940 6002 94 31533 31627 
31 869 900 63 5940 6003 95 |>63324 |>63419 




















32 868 900 64 | 5940 6004 | 96 — 
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Calculating the Coefficients of Certain Linear 


Predictors 


By D. S. Stoller and L. C. Stoller 


It is assumed that observations, x; , have been made at the n + 1 points, 7 = 


0, 1, --- , n, which are equally spaced. It is desired to find a linear predictor 


(1) Yn41 = AX + ore + Ann 





Received 28 September 1958. 
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such that yn4: is an estimate of the “next” observation, z,,, , which has not yet 
been made, under the assumptions* 


(a) om =***-=o,,=0 


(2) (b) E(2z;) is a k* degree polynomial in j 
(ec) E(yas) = E(2n4:) for an arbitrary polynomial of degree k. 


As an ancillary result, Karush and Wolfsohn [1] obtain expressions for the coeffi- 
cients, a;, which minimize o;,,,/0°. These expressions involve determinants of 
order k. Contained herein is an expression for these coefficients which is more 
convenient for computational purposes. Also a table of some of these a; is deposited 
in the unpublished mathematical tables repository. (See Reviews and Descriptions 
of Tables and Books.) By assumption (b) in (2), one can write 


(3) E(x;) = aoPo(j) + mPi(j) + --- + oePe(7); jG =0,1,---,n +1 


where the a; are constants and where the P;(j) are polynomials in j of degree i 


which are orthogonal over the set of points, (0, 1, --- , 7), Le., 
(4) 2d P.(s)Pc(j) = 0 st. 


Taking expected values on both sides of equation (1), one can write 
(5) E(Xn41) = oH (xo) + +--+ + Gn (tn) 
or, from (3), 
aoPo(n + 1) +--+ + oxPi(n + 1) 
= aofaoPo(0) + --- + axP(0)] + --- + aalaoPo(n) + --- + oxPr(n)]. 


From (6) one obtains the following set of equations which the coefficients must 
satisfy : 
aPo(0) + --- + anPo(n) = Po(n + 1) 


aoP.(0) + --- + anPi(n) = Pi(n + 1). 
Since q 
(8) “ust = ag’ + +++ tae, 
it is desired to minimize (8) under the constraints (7). Form the expression 
F = (ag + +++ + ax) + dofaoPo(0) + --- + anPo(n)] 
+ +++ + AfaoP.(0) + --- + anPi(n)] 
where the \; are Lagrange multipliers. Setting dF /da; = 0, one obtains 
(10) 2a; + AoPo(j) + +++ + AUP. (7) = 0 j = 0,1,---,n. 


*The symbol, oz;, refers to the standard deviation of z;; ¢ is a constant; and E(z;) refers 
to the expected value of 2; . 


(9) 
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If these equations are all multiplied by Po(j) and then summed, and this process 
is repeated but with P;(j), --- , Pi(j), one obtains the equations (by using (4) 
and (7)): 7 
2Po(n +1) + rod Pe(j) = 
. ~— 
(11) 
2Pi(n +1) + me Do Pi'(j) = 0. : 
i= 
Therefore, by substituting for the A; in (10), one obtains 
m Po(j)Poln + 1) eh + Pl Paln + 1). 
dX Po'(3) D Pi'(j) 








(12) = 0,1, --- ,n. 


7=0 7=0 ] 


Values for a; were obtained by using the following set of polynomials, which are 
orthogonal over any set of points symmetric about the origin [2]: 


Py(t) = 
P(t) =¢ 
(13) 
Pi4,(t) = tP;(t) = BiPia(t) 
where . j 
(14) Bi = > P?(t;) > Pi_-1(t;). : 


By transforming the integer set 0, 1, --- , m to the interval (—2, 2) we derive [3] 
a simple expression for the (; : 


4i'[(n + 1)? — 7] 
n?(47? — 1) : 

It may be noted that >-}-. a; should be equal to one. For degrees one through 

eight this check sum for the coefficients is verified with a minimum of six signifi- 


cant figures. For degrees nine and ten the check sum is verified to a minimum of 
five significant figures. 


(15) Bi = 





$=1,2,---,k—1. 


RAND Corporation, 
Santa Monica, California and 
Space Technolo Laboratories, 





Los Angeles, California 
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REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


13([D, Z].—E. G. KoGBer.iantz, “Computation of Aresin N for 0 < N < 1 using 
an electronic computer,” IBM, Journal of Research, v. 2, No. 3, 1958, p. 218-222. 


All known subroutines for Aresine are based on the relation Arcsin N = Arctan 
N/1i-wN ie Therefore, Arcsine is not computed as such, but as an Arctangent. 

To avoid the loss of machine time caused by the computation of N/(1 — N*)', 
a direct computation of Arcsine is proposed. A subroutine yielding the first six correct 
significant digits in only five multiplications and divisions is described in detail to 
illustrate the new method’s rapidity. The same number of five operations is neces- 
sary to compute, knowing N, the number N/(1 — N’)?. 

AUTHOR’s SUMMARY 


14[D, Z].—E. G. KoGBer.iantz, “Computation of Arctan N for —«» <N < + 
using an electronic computer,” IBM, Journal of Research, v. 2, No. 1, 1958, p. 
43-53. 

Rational (R) and polynomial (P) approximations to Arctan N are studied with 
the aim of computing this function, to any prescribed accuracy and without unduly 
increasing the number PC of stored constants, in a minimum number M of multi- 
plications (and divisions for R approximations). The number Dg of first correct 
significant digits in principle is not bounded. The results corresponding to the values 
8, 10, 18 and 20 of this number are as follows: 











Single Precision Double Precision 
Approximation Point (Computation) 
Dg M PC Deg M PC 
(4 21 (6 19 
Floating 8 } 18 
5 9 \7 17 
Rational (R) 
| (5 14 (6 30 
Fixed 10 20 { 
6 9 [7 18 
(5 10 17 en 
Floating 8 
\6 8 18 9 14 
Polynomial (P) 
| (6 1 [eo @ 
lFixed 10 20 
\7 9 10 15 


If M is increased, subroutines with smaller PC are easily deduced from our 
general results. Thus, for instance, rational approximations with Dg = 6 can be ob- 
tained in three multiplications only, if PC = 19, but the same accuracy Dg = 6 
characterizes also the cases M = 4 with PC = 11, and M = 5 with PC = 7. 

If polynomial approximations are used, Dg = 6 is obtained for M = 5, PC = 7, 
but also for M = 4 and PC = 11. No subroutines with a stored table of values of 
Arctan x are considered. 

AvuTHOR’s SUMMARY 
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15{E, Z].—E. G. Koeperuiantz, “Computation of e” for —0o < N < + using 
an electronic computer,” IBM, Journal of Research, v. 1, No. 2, 1957, p. 110- 
115. 


Rational (R) and polynomial (P) approximations to the exponential function 
e” are studied. They allow e” to be computed for any value of the exponent N in 
the infinite range from minus infinity to plus infinity in a minimum number M of 
multiplications (and divisions, for the rational approximation). This minimum is 
attained without unduly increasing the number PC of precomputed and stored 
constants and also without limiting the number Dg of the first correct significant 
digits. The main results are presented in the following table. 





Single Precision Double Precision 
Machine Approx. Computation = ~ De 7 ~; D ~ 
ee point 3 19 9 5 21 19 
Binary R 4 8 10 6 8 18 
\fixed point 3 35 10 6 10 21 
4 8 10 
(floating point 5 17 8 8 30 18 
Decimal a 9 24 17 
\fixed point 5 25 10 9 32 20 
6 19 10 


AvuTHOR’s SUMMARY 


16[F]—A. GuopENn, Table de Factorisation des Nombres N* + 1 dans l’intervalle 
10000 < N S 20000. Manuscript of 50 leaves deposited in UMT File. 


The present table completes that of the same title reviewed in MTAC, v. IX, 
p. 26. In the present one, thanks to the collaboration of Dr. N. G. W. H. Beeger, 
a considerable number of large factors less than 64-10", identified as prime num- 
bers, have been entered. 

The machine calculation was performed in the following manner. The number 
N* + 1 or (N* + 1)/2 was divided by the product of the known prime factors, 
then the quotient was divided successively by each of these factors to ascertain if 
any was a multiple factor. If none of these last divisions was exact, the quotient 
was identified as a prime number. 

Some typographical errors appearing in the previous manuscript table have been 
eliminated. 

My “Table des solution de la congruence z* + 1 = 0, mod p, pour 8-10° < 
p < 10°” is in preparation. This will permit the entry in the present table of the 
prime factors occurring in that interval. 

A. GLODEN 


11 rue Jean Jaures 
Luxembourg 


17(F).—Hansras Gupta, C. E. Gwytrner & J. C. P. Mituer, Tables of parti- 
tions, Royal Society Mathematical Tables, v. 4, Cambridge University Press, 
1958, xxxix + 132 p., 28 em. Price $12.50. 
The partition function p(n, m) = p,(n, m) is defined as the number of partitions 
of n into parts not exceeding m. More generally p,(n, m) for s > 1 is the number of 
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partitions of n into parts where there are s varieties of each of the integers 1, --- , m 
and s — 1 varieties of the integers greater than m. 

These tables are four in number, giving respectively the values of p(n, m), 
p(n, m), ps(n, m) and p,(n, 0). Each table provides a partial check on the preced- 
ing table. Thus, the formula 


P(n,m) = > (— 1)’pe (n - (m+ (E>) + 3) +7) 


where R is the largest value of r for which the first parameter is not negative, pro- 
vides a check on Table I from Table II, and may be used to calculate values of 
p(n, m) outside the range of Table I. 

The tables are preceded by a fairly long introduction and an extensive bibliog- 
raphy. The introduction describes the tables, gives historical background, the basic 
recursions, and exact formulae for p(n, m), m = 1, --- , 12. The exact formula for 
p(n, m) is a polynomial of degree m — 1 in n, plus terms of lower degree depending 
on the residues of n with respect to moduli not exceeding m. The exact formulae 
involve inconveniently large numbers as is illustrated by the calculation of p(100, 
12) = 167 13148. Theoretical discussions and some asymptotic formulae are 
included. 

Table I gives p(n, m) for n from 1 to 200 with m ranging from 0 to min (n, 100), 
and from n = 201 to 400 with m ranging from 0 to 50. Table II gives p.(n, m) for 
n from 1 to 1000 with m ranging from 0 to the value indicated by max m below: 


Range of n 1-50 50-100 100-150 

max m n 23 20 

Range of n 150-200 200-250 250-300 

max m 12 11 7 

Range of n 300-350 350—400 400-450 450-500 
max ™ 6 5 4 3 
Range of n 500-550 550-1000 

max m 1 0 


Table III gives p3(n, m) 


for n = 1 to 50 withO <mazn 

forn = 50to 100 withO < m < 19 

forn = 100 to 150 withO << m <6 
and forn = 150 to 200 with m = 0. 


Table IV gives p(n, 0) for n = 1 to 200. 


MarsHALL Hatt, Jr. 
Ohio State University 
Columbus, Ohio 


18{I, X|—Her Masesty’s Nautica AtMANnac Orrice, Subtabulation, A Com- 
panion Booklet to Interpolation and Allied Tables, Her Majesty’s Stationery 
Office, London, England, 1958, 54 p., 24.5 em. Price 7s.6d. net. 


Interpolation and Allied Tables was reviewed recently in this journal [1], and it 
was then pointed out that it contained material essential for every computer and 
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numerical analyst, whether practicing or aspiring. This booklet, in which the print- 
ing and layout are excellent, contains material concerning methods of systematically 
interpolating a function, tabulated at an interval h, tc obtain its values at a smaller 
interval, h/p, where p is an integer. Since this is a problem which is almost certain 
to be encountered by computers early in their career, they should be aware of 
efficient ways of dealing with it. Three types of methods are described, and their 
merits under various conditions, principally with reference to equipment available, 
are discussed. 


I. Direct Methods 


These methods are the most straightforward: one simply interpolates using his 
favorite method (or the one most suitable in the circumstances). We describe the 
tables which are provided, using the notation of the booklets, which is summarized 
in the earlier review [1]. 

1. Bessel’s Formula 

The coefficients to 6D of the double second difference, and of the third and fourth 
difference corrections are given in the cases p = 20, p = 24. [These cases, of course, 
cover those when p is a factor of 20, or 24.] 

2. Lagrange’s Formula 

Six-point coefficients, to 8D, are given in the cases of p = 20, p = 24. 

3. Everett’s Formula 

Again for p = 20, p = 24, the coefficients Z, and F2 are given to 7D, M, and 
N, to 3D, E, and F, to 6D. 

4. Lagrange’s Formula 

Ten-point coefficients, to 10D, are given in the cases p = 20, p = 24. 

5. High-precision formulae and coefficients 

These cover the cases p = 2, p = 3 and p = 10. 

6. Throw-back formulae 

This is a useful collection of formulae covering various simple, simultaneous 
and generalized throw backs for the Bessel and Everett formulae. 

There are various worked examples and sound, practical advice on detailed 
points. 


II. Precalculated second differences 


This method is a new development which supersedes the ‘‘end-figure” method 
of Comrie, and is suitable when only desk machines are available and in cases 
where the sixth difference can be as much as 500,000. 

The basic idea is to take the Bessel interpolation in the form 


fo = fo + p83 + Bo(Smo + dma) + Bodm + Ta(yo' + v1) + Tory 
where, in addition to the notation of the previous review, we have used 
1000y° = 6 — 0.26’ + --- 
Ts = 10000(B; + 0.108B;). 


This can be written as 
fo = {fo + rA + r(r — n)B} + {4vr(r — n)(2r — n)C} 
+ {rva + Bo(rv)b} + By(rv)C + Tayo? + 11’) + Tor 
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where A, B, C, a, b, c are integers defined by 
4=nA+a, (in +im) =4n°B+b, bm = 2n'C +c, 


and where p = r/n = rv, r being an integer. 

Then the terms in the first braces have a constant second difference (with re- 
spect to r) of 2B; the term in the second braces has a constant third difference of 
2C. Tables are given to provide the second differences of the remaining terms, and 
the leading first differences. From these the values of f, can be built up. There are 
worked examples, illustrating two methods of use of the tables. 

The tables provided cover the cases p = 4, 5, 6, 10. 

These methods, although very powerful, are rather complicated, and are not 
likely to appeal to the more casual user, who will probably be content with the 
methods described in the first section. 


III. The method of bridging differences 


This is a method which is appropriate for use with multi-register adding ma- 
chines. There is a discussion of the theory, and a list of formula indicating their 
limits of applicability and the number of registers needed. 


Joun Topp 
California Institute of Technology 
Pasadena, California. 
1. MTAC, Rev. 54, v. 12, 1958, p.99-103. 


19/K].—D. 8S. Srotter & L. C. Stouuer, Tables of the Coefficients of Certain Linear 


Predictors. 
Tables are given for the coefficients ap , a; , --- @» of the linear predictor yn4; = 
Q Xo + +++ + Gntn, Where 2%, 21, --* Z, are Observations atn + 1 equally spaced 


points and y,4: is an estimate of the next observation, z,4; . Tables are listed for 
polynomials of degree 1 to 10 inj of the expected value H(2z;) of x; and for var- 
ious numbers of observed points. A detailed description of the assumptions under 
which these tables are computed and discussion of the method of computation are 
given in [1]. These tables are on file in the Unpublished Mathematical Tables reposi- 
tory of MTAC. 

7. 


1. D. 8. Srotter & L. C. Srouusr, “Calculating the Coefficients of Certain Linear Pre- 
dictors,’”” MTAC, v. XIII, number 66, April 1959, p. 122-124. 


20{L].—Carson Fitammer, Spheroidal Wave Functions, Stanford University Press, 
Stanford, Calif., 1957, ix + 220 p., 25.5 em. Price $8.50. 


A condensed, but fairly complete collection of definitions of the angle and radial 
solutions of the wave equation in prolate and oblate spheroidal coordinates is pre- 
sented in this book, together with a set of tables relating the notations for the 
various related functions, which are now in use in the literature. Series expansions 
and various integral representations for the functions are given. A bibliography of 
work on these functions is given, and the more important results are brought into 
line with the definitions used in the present book. A useful chapter is included on 
the related solutions of the vector wave equation in spheroidal coordinates. 
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The numerical tables have been adopted from several sources, reworked, and 
previous errors rectified. Coefficients in the series of powers of c (interfocal distance 
divided by twice the wave length) for the separation constants A, and for the co- 
efficients d-” of the expansion of the prolate angle functions in associated Legendre 
functions P®+,-(7) are given to ten decimal places for m = 0(1)3, n = m(1)9. The 
separation constants themselves are given to at least six decimals for m = 0(1)3, 
n = m(1)3, c = 0(.2)5, and for m = 1, m = 5(2)19, for a few values of c between 
1.2 and 3.2, including $7, $7 and z. A similar tabulation is given of the expansion co- 
efficients d for the angle and radial functions, and also for the coefficients of the 
expansion of the angle functions in powers of (1 — 7’) (» being the angle coordinate 
in the elliptic coordinates which form the spheroids by axial revolution). 

The prolate angle functions are also tabulated, to four decimals, for m = 0(1)3, 
n = m(1)3 for c = 0(.5)5 and 0(= cos‘) = 0(5°)90°. Values of the prolate 
radial functions of the first and second kind, and their derivatives are also given, 
for these values of m, n and c, for a few values of the radial coordinate ~ near the 
origin ( = 1). 

Corresponding values for the oblate functions, for a somewhat reduced range of 
values of c, m and n are also tabulated. 

The book should be welcomed by those working in the field. However, this 
reviewer is sure the author will agree that more work, both analytic and computa- 
tional, is needed before wave calculations in spheroidal coordinates will be com- 
paratively easy to carry out. 

Puitie M. Morse 


Massachusetts Institute of Technology 
Cambridge, Massachusetts 


21(L].—M. E. Lynam, Table of Legendre Functions for Complex Arguments, The 
Johns Hopkins University Applied Physics Laboratory Report, TG-323, 1958, 
4 p., Tables. 


The Legendre functions P,(x) and Q,(x) and their first derivatives are tabu- 
lated to either 8D or 8S for z = 0 (.05) .95 and for values of » with negative im- 
aginary part that satisfy the equation v (1 + v) = —zy, where y = 0.01 (.01) .10, 
0.2 (.1) 2.5, and 3 (1) 24. 

An indication of the relative accuracy of the tabulated values is given in a 
double column headed “Error,”’ which records the (negative) characteristics of the 
common logarithms of both the real and the imaginary parts of the difference 


1-(i- | P(e) f Q(z) — Q(x) 2 pz) | 


when evaluated by substituting the computed values in the Wronskian determinant. 
Spurious symbols in this column are due to printing difficulties arising from char- 
acteristics numerically greater than 9. 

The author states in the text that the calculations were performed on the Univac 
Scientific 1103A Computer, employing a subroutine previously prepared for calcu- 
lations of rocket instability. 


J. W. W. 





Sp PRAY ae 


ETE NS OLS 








22(L 


pure 
Sn(2 
ing, 


Scie 


give 
sing] 
On | 
valu 


pole 
whic 


mod 


to ; 
tabl 





._ oo. ie J SE eed 


£ 


Ow ODO ww 


f 


\- 
\- 





7 TENE Se 








REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 131 


22{L].—Morio Onog, Tables of Modified Quotients of Bessel Functions of the First 
Kind for Real and Imaginary Arguments, Columbia University Press, New York, 
1958, ix + 338 p., 27.3 em. Price $12.50. 


As is generally known, tables of Bessel functions of integral order for real and 
pure imaginary arguments are now adequate. But since the “modified quotient” 
,(z) = zJn.+(z)/J,(z) frequently occurs in mathematical physics and engineer- 
ing, this book gives tables of $,.(z) and $, (ix), in both cases for n = 1(1)16 and 
x = 0(.01)20, without differences. Most of the computation was done at the Watson 
Scientific Computation Laboratory, chiefly on an IBM 607. 

The tables of %,,(x) occupy pages 17-176. The number of significant figures 
given is normally between eight and ten, and falls below six only in a few cases at 
single arguments near zeros (or, in one instance, a pole) of the functions tabulated. 
On page 174, it may save trouble to note that at n = 14, z = 18.90, the printed 
value 92338, followed by two spaces and a decimal point, means 9.2338- 10°. 

The tables of %,(ix) occupy pages 179-338. As the functions have no zeros or 
poles in the range of tabulation, the values are given uniformly to eight decimals, 
which means that there are either nine or ten significant figures. 

The short introductory text includes a number of formulas relating to the 
modified quotient functions. 

A. F. 


23{L].—James C. Wittse & Marcia J. Kine, Values of the Mathieu Functions, 
The Johns Hopkins University Radiation Laboratory, Tech. Report No. AF-53, 
Baltimore, Maryland, 1958, 78 p. 


{L].—James C. Wuatse & Marcia J. Kina, Derivatives, Zeros, and Other Data 
Pertaining to Mathieu Functions, The Johns Hopkins University Radiation 
Laboratory, Tech. Report No. AF-57, Baltimore, Maryland, 1958, 75 p. 


The authors have used the notation and normalizations defined in [1]. In the 
following summary, notation such as r = 0(1)2 will imply that r = 0 does not apply 
to functions associated with odd periodic solutions, such as So,(s, x), Jo,(s, x). 

The following tables appear in the first report under review: 

Table I. Periodic Mathieu functions Se,(s, v) and So,(s, v), r = 0(1)2; 
s = 3, 6, 9, 15, 25, 40; 7 to 13 angles v between 0° and 90°; for Se.(s, v), the addi- 
tional parameters s = 1(1)9; Se,(—s, v) and So,(—s, v), r = 0,1; s = 3, 6, 9, 15, 
v = 0(15°)90°; 4D. 

Table II. Radial Mathieu functions of the first kind, Je,(s, uw) and Jo,(s, u): 
r = 0(1)2; s = 1, 3, 6, 9, 15, 25, 40, for 15 to 22 values of u less than 2.01, un- 
equal intervals in u; 3D or 4D. 

Table III. Radial Mathieu functions of second kind, Ne,(s, u), No,(s, uw), s up 
to and including 25; same range of u and format as Table IT. 

Table IV. Mathieu-Hankel functions He,”’(—s, u) and Ho,” (—s, u), s = 1, 
3, 6, 9, 15; r = 0(1)2; 2D to 4D. Format and range of u about as in Table II: for 
r = Oand 2, s = 25 also. 

(In all tabulations in both reports, values of the function at s = 0 are given 
to 5 or more significant figures, since they could be readily obtained from the 
tables in [1]. ) 
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An important feature of the report is a set of 39 pages of good graphs from which 
one may in fact read the magnitude of the functions or better, in the region covered 
by the tabulation. Quoting the authors: “The tables contain approximately 1750 
calculated values of the various Mathieu functions.” 

The following tables appear in the second report: 

Table I. Se,’(s, v) and So,’(s, v), r = O(1)2, s = 1, 3, 6, 9, 15, 25, 40; » = 
0(15°)90°; 4D. Also Se,’(—s, v), So,’(—s, v), s = 1, 3, 6, 9, 15; »v = 0(15°)90°; 4D. 

Table II. Je,’(s, u) and Jo,’(s, u); r = 0(1)2; s = 1, 3, 6, 9, 15; for r = 1, 
s = 25, 40 also; 2D-4D. Similar range of u to that in Table II of the first report. 

Table III. Ne,’(s, u), No,’(s, u), r = 0(1)2, s up to and including 15; 2D to 
3D. 

Table IV. Derivatives with respect to u of He,” (—s, u) and Ho," (—s, u); 
2D-4D; same format and range of u as Table III. 

Table V. Seo(—s, v) and So.(—s, v); s = 1, 3, 6, 9, 15; v = 0(15°)90°; 4D. 

Table VI. First Zero of Je,(s, u) and Jo,(s, uw), other than u = 0;r = 0, 1, 2; 
s = 3, 6, 9, 15, 25, 40; also s = 1 for r = 0; 3D. 

Table VII. Second Zero of Je,(s, uw) and Jo,(s, u); s = 6, 9, 15, 25, 40; for r = 
0, also s = 3; 3D. 

Table VIII. First Zero of Je,’(s, u) and of Jo,’(s, u), r = 0, 1, 2; s = 3, 6, 9, 
15; for r = 1, also for s = 1, 25, 40; for r = 2, also s = 1; 3D. 

The tables are followed by 32 pages of graphs of the various functions. 

The periodic Mathieu functions for the parameter s can be found in the very 
accurate and systematic tabulation of Ince [2]. However, the authors devoted rela- 
tively little space to them, and since the normalization and notation are different, 
there was every justification for including this material. The tables for the param- 
eters —s and the radial functions are new, and form a worthwhile exploratory 
contribution. The graphs in both pamphlets constitute an especially welcome 
addition. 

The writer spot-checked two values of the radial functions from some systematic 
tabulations that are being made at WADC. The values are essentially as accurate 
as claimed, although the authors’ argument, such as u = 0.182, does not mean that 
the authors used precisely this value of u, but one close to it, which gave them a 
convenient value for subsidiary arguments required in their computations. It would 
be too difficult to check the accuracy of the tables systematically, since the u-argu- 
ments are not equally spaced. Some values of Ne, (s, u) on page 29 of the first report 
have been changed by pen in the reviewer’s copy. However, the tables were not 
intended for interpolation purposes, and even if some other errors still exist there, 
the tables serve a useful purpose for exploratory studies. The reviewer noted only 
a few inaccuracies: 

On pages 6 and 7 of the first report, the normalizations as given are not com- 
pletely defined. A better statement would have been, “... the functions are nor- 
malized so as to satisfy equations (12)—(17).” 

In Table IV of the second report, corresponding to s = 9, the “‘i’’ was omitted 
on the first line. 

Two sentences on page 9, second paragraph of this report need revision. The 
radial Mathieu functions approach the corresponding Bessel functions, as s ap- 
proaches zero, only if u approaches infinity in a special manner. It is not true that 
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the Mathieu functions approach the Bessel functions, as s approaches zero, when 
u is merely different from zero. The special manner requires that s’ cosh u approach 
a constant value kp, in the authors’ notation. A similar change is necessary on page 
9 of the first report, first paragraph. 

The authors deserve credit for the careful performance of a very useful task that 
presented considerable difficulty. 


GERTRUDE BLANCH 
Aeronautical Research Center 
Wright Air Development Center 
Wright-Patterson Air Force Base, Ohio 


1. NBS, Tables Relating to Mathieu Functions, Columbia University Press, New York, 1951. 
2. E. L. Ince, ‘“‘Tables of Elliptic Cylinder Functions,”’ Roy. Soc. Edin., Proc. 52, 1932, 
p. 355-423. 


24[L,S|—-W. M. Rocers & R. L. Poweut, Tables of Transport Integrals 
Jn(x) = eae NBS Circular 595, Government Printing Office, Washing- 
5 ~ 
ton, D. C., 1958, ii + 46 p., 26 em. Price 40 cents. 


The transport integrals J,,{«) defined in the title may be expressed in terms of 
the integrals 
ae 
e—l 





F,(x) = I 


by the formula J,(z) = nF,4(x) — x”/(e* — 1). Tables of one kind or another 
have been made for various values of n by various authors, but a consolidated table 
has been badly needed and the appearance of one is welcome. 

The main tables (pages 7-45) give values of J,,(z) to 6S without differences for 
n = 2(1)17, x = 0(.1)X, where X depends on n in the following way: 


n 2(1)5 6, 7 8(1)11 12(1)17 
x 25 30 35 40 


For n = 2(1)7, 6-figure values of J,(x)/x"” are also given. 

The values of J,( ©) may be expressed in terms of the Riemann zeta function 
by the formula J,( ©) = n!¢(n). Small tables on page 46 include values of ¢(n) 
and J,() to 8S for nm = 2(1)17, as well as values of Bernoulli numbers used in 
expansions. 

The computations were performed on an IBM 650. Details are given regarding 
accuracy, which varies with n and x. Usually the maximum error appears to be 
either one or two units of the sixth figure, but it seems that it may be as much as 
five units for nm = 12(1)17 in the range of z from 9 to 11. Proper understanding is, 
however, made difficult by instability of notation, in that series expansions for small, 
intermediate, and large values of x are respectively referred to sometimes as a, b, c, 
(for example, in Table 1) and sometimes as c, b, a (for example, on pages 2-3). 

The introduction gives a mathematical account of the functions tabulated, along 
with references to earlier tables. To these it seems worth while to add the tables of 
Nagai and Umeda [1], which, as far as J,(2z) is concerned, give 6S for n = 5(2)11, 
xz = 0(.1)10. These tables are in internationally intelligible form, though the ac- 
companying text is in Japanese. Comparison between the two tables shows that 
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Nagai and Umeda’s value at n = 9, x = 0.8 is entirely wrong. Apart from this, 
differences rarely exceed two units of the sixth significant figure, and never exceed 
five units. Such of the larger discrepancies as the reviewer investigated all turned 
out to be due to errors in the Japanese tables, which date from the pre-automatic age. 


A. F. 


1.8. Nagar & K. Umepa, ‘‘Rikwagaku-kenkyi-jo ihé,’’ Bulletin of the Institute of Physical 
and Chemical Research, v. 18, no. 7, 1939, p. 529-536. 


25(M, Z].—A. R. DiDonato & A. V. Hersuey, New Formulae for Computing In- 
complete Elliptic Integrals of the First and Second Kind, NPG Report No. 1618, 
NAVORD Report No. 5906, January, 1959. 


The authors study the application of series expansions for efficient evaluation of 
elliptic integrals of the first and second kind on high speed computers, with special 
reference to the NORC. To facilitate computation when the modulus is near unity, 
they derive an appropriate series expansion for each incomplete integral. The for- 
mulae are essentially of the same character as previously given by Radon [1], who 
has made a rather thorough investigation of the subject. Radon gives numerous 
expansions and also studies the incomplete elliptic integral or the third kind. 


YupeEti L. LuKE 
Midwest Research Institute 
Kansas City, Missouri 


1. Briertre Rapon, “Sviluppi in serie degli intergrali ellittici,’? Accad. Naz. Lincei, 
Atti Mem., Cl. Sci. Fis. Mat., s. 8, v. 2, 1950, p. 69-109. See also MTAC, v. V, 1951, p. 79-80. 


26(P].—Rene A. Hiconnet & RENE A. Grea, Logical Design of Electrical Circuits, 

McGraw-Hill Book Company, Inc., New York, 1958, ix + 220 p., 24 cm. Price 

$10.00. 

This book is devoted to a binary (Boolean) analysis of the synthesis of switching 
circuits. Most attention is given to relays and relay circuits, but there is a chapter on 
vacuum-tube and diode circuits. In the main, the book is a convincing and realistic 
discussion by people experienced in relay circuitry. 

The notations explained and used include Euler circles, ordinary Boolean algebra 
(using the plus symbol + for ‘‘or’, either a multiplication dot or more frequently 
no mark for “and,” and a prime ’ for “‘not’’), and reference to the vertices, edges, 
and other parts of the boundary of the unit cube in n-dimensional space. It seems 
unfortunate that books dealing with this subject continue to use the Boolean sym- 
bols which can be confused with ordinary arithmetic symbols instead of v and A 
for “or” and “and”’ respectively, but this continues to be common. 

Many examples make the book easy to follow. However, absence of statements 
of theorems or other displayed punch lines will force most readers at least to skim 
through the entire book rather than to refer to isolated sections in order to get 
much from it. 

An appendix includes a table of four-relay contact networks prepared by Edward 
F. Moore. 
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There are no bibliographic references not immediately incident to the material 
presented. 
The book is a revision and translation of an earlier French edition. 


a 


C. B. T. 
27(S).—James RippE.1, A Table of Levy’s Empirical Atomic Masses, Report AECL 
i No. 339, 195€, 89 p., 27 em. Available from Scientific Document Distribution 


, Office, Atomic Energy of Canada Limited, Chalk River, Ontario, Canada. Price 
7 $2.00. 
Nearly 4,000 atomic masses have been calculated from Levy’s empirical mass 
| equation. These, together with the neutron, proton, and alpha-particle binding 
» & energies which can be calculated from them, are tabulated in the report. 


AvuTHOR’s SUMMARY 


28[W, Z].—DrEPaRTMENT OF THE ARMY PAMPHLET, Introduction to Automatic Data 
Processing, No. 1-250-3, Headquarters, Department of the Army, Washington 
25, D. C., 1958, 88 p. 


This pamphlet is an introduction to the concepts, methods, and equipment 
components of automatic data processing systems applied to management problems. 
A large part of the material is devoted to applications with emphasis on the utiliza- 
tion of automatic data processing facilities by the U. S. Army. 

Various types of data processing equipment are introduced in the first chapter. 


= 
nm 
TTL SS A SEMA on ke 





No previous knowledge of computers on the part of the reader is assumed. After an 
i, FF explanation of the concepts of machine language and the binary number system, a 
P meaningful description of the more common storage devices follows. The pertinent 
E physical characteristics of the different components are concisely stated with in- 
Is, formative comments on the advantages and disadvantages of each. This section is 
” replete with illustrations and diagrams of typical computer systems and compo- 
nents. 
8 The fundamentals of programming are also included in the comprehensive first 
- chapter. This presentation is by no means a complete description of the discipline 
- of computer programming, but it is a straightforward explanation in simple terms 
of the most important principles. 
™ The title of the second chapter is “Management Aspects of Automatic Data 
ly Processing Systems’’. The chief concern of management is neatly stated as the “con- 
> trol of men, money, minutes, materials, and machines”. The thesis here is that all 
” this can be done most effectively by the intelligent use of automatic data processing 
wi systems. To guide managers who are about to consider such a system, timely advice 
. is given on cost considerations and on the difficulties of evaluating the production 
gains to be realized once the system is installed. The questions of purchase versus 
on rental, site planning, and personnel requirements are treated realistically. 
= The problem of communication with the heart of a data processing system (i.e., 
et with the main computer) is considered in the third chapter. Whenever information 
“a is to be transmitted to and from widely dispersed operation centers, fast and reliable 


communication devices must be selected for use in the system. Data transmission 
devices in current use are described briefly, and examples of how they are incorpo- 
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rated into large systems are to be found in the fourth chapter, entitled “Applica- 
tions of Automatic Data Processing Systems’’. The varied applications described 
are examples of actual systems in operation taken from private industry, govern- 
ment, and the U. 8. Army. 

The fifth and final chapter elaborates further on the Army utilization of auto- 
matic data processing systems. A number of proposed uses indicate the trend the 
Army is likely to follow in extending automation to its many activities. 

This booklet would be more valuable if it included a bibliography of its source 
material as well as suggested references for additional information, particularly on 
equipment and programming. Nevertheless, it is a useful presentation of the actual 
and potential benefits to be gained, as well as the problems to be faced, in employing 
automatic data processing systems in large-scale organizations whose operations are 
widespread both geographically and functionally. 


J. W. Scuor 
Applied Mathematics Laboratory, 
David Taylor Model Basin, 
Washington, District of Columbia 


29(X, Z)|.—K. H. Booru, Programming for an Automatic Digital Calculator, Aca- 
demic Press Inc., New York, 1958, vii + 238 p., 22 em. Price $7.50. 

The value of this book lies in its provision of many examples of codes for com- 
mon operations using the automatic computer APEXC. Generally, these operations 
are applicable to other computers. There are several major weaknesses, however, 
from the point of view of the reviewer. Most important is lack of precision of lan- 
guage and symbolism. Next in importance is lack of an adequate bibliography or 
reference to methods other than those recommended for the particular machine; 
this paucity of references has an irritating impact on the reviewer also in that he 
finds no credit to many of the early workers in the field. A third weakness lies in the 
completely pervading reference to the APEXC machine, which is not typical of 
most of the machines presently available throughout the world. The book is a good 
coding manual for the APEXC machine, but it falls short of being a well rounded 
work on coding and its logic. 

Chapter 1 is introductory in nature, and describes the machine. Chapter 2 
describes the techniques of programming (our usual terms of programming and 
coding, flow chart, and so on are not used, but the author uses the ‘schematic pro- 
gramme” and the “detailed programme”’). Here there is some material that is ap- 
plicable only to machines with cyclic memory reference and with provision for op- 
timal spacing of commands in the sequence. 

Chapter 3 discusses input and output, and again it is highly directed toward the 
one machine. It contains a table of decimal-binary equivalents for integers n = 
0(1)31 (decimal). There is a detailed program for input of binary data, an input 
check program, a conversion program for decimal numbers, a decimal output pro- 
gram, a binary output program, and a program for storing data from decimal input. 

Chapter 4 is devoted to division and square root. A division routine which de- 
velops the quotient in digits —1 and 1 instead of 0 and 1 (nonrestoring division) 
and then corrects is included. The usual iteration « — (x + b/x)/2 for the square 
root of b and the usual iteration for the inverse square root (including normaliza- 
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tion) are coded. Also coded are a method of successive subtraction of odd numbers 
for square roots, and a bisection method for the extraction of higher roots. 

Chapter 5 includes codes for the evaluation of polynomials (using Horner’s 
method), of even polynomials, and of odd polynomials; it contains explicit codes 
for the evaluation of the sine and cosine functions by Hastings’ polynomial ap- 
proximation [1], a code for tabulating a function with constant n-th differences, and 
a code for Simpson integration using no more than 31 ordinates. 

Chapter 6 contains some “‘non-arithmetic codes’’. These include codes for nor- 
malizing numbers, selection of the largest modulus from a set of numbers, and col- 
lating. 

Chapter 7 is devoted to matrix operations. It includes codes for matrix multipli- 
cation and for finding the largest real eigenvalue of a matrix by the power method. 

Chapter 8 has a code for solving no more than 20 linear algebraic equations in 
20 unknowns simultaneously, using the elimination method. 

Chapter 9 is devoted to machine translation from French to English, using stored 
stems and suffixes. Due account is accorded to the possibility that the adjective 
following the noun in French should precede it in English. 

Chapter 10 returns to the interpretive program for decimal tapes. This gives in 
some detail one of the conversion difficulties which beset the user of a binary ma- 
chine. 

Chapter 11 is devoted to interpretive routines and pseudo-codes. Floating-point 
arithmetic is the only subject treated in detail, and codes for the elementary opera- 
tions are included. 

Chapter 12 is devoted to checking. It includes ordinary post-mortem routines. 

An appendix includes some well selected minor tables. Here one finds the first 
thirty powers of two expressed decimally, the first nine positive and negative powers 
of ten expressed in octal digits, the binary equivalents of the simplest decimal frac- 
tions, and octal expressions for one hundred twenty-six 18 or 2D decimal numbers. 

Thus we have an instructive collection of codes available, and this collection 
should be valuable to anyone working with a new computer, anyone who has access 
to a computing center with similar codes not carefully explained, or anyone who 
wishes to see the current practice of explaining and printing codes. In general, this 
presentation is good. 

Now, without belaboring the point, some of the less agreeable features of the 
book should be mentioned. 

The book is really more a collection of codes than a book on programming, so it 
may be misnamed. No mention is made of the underlying logical principles of cod- 
ing; however, these have not been presented successfully anywhere else. Still, the 
inductive process in coding (or the use of difference equations or recursive relations) 
is not stressed at all. 

There is minimal mention of the work of others. This is important in that the 
reader is not given clues as to other places to look for information. There seems to 
be reason to believe that the author is unaware of much pertinent literature; her 
complete ignoring of the work of others in machine translation makes this chapter 
almost useless, even though (as she admits) her group is a leading center in this 
work. 
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The author’s language and notation are imprecise, and it seems that the neo- 
phyte coder might deduce that precision is not necessary (a position which seems 
completely untenable, since coding is nothing more than the precise expression of 
arithmetic procedures in a restricted basic language). For example, on page 113 
there is a careless confusion in subscripts; in Chapter 5 it is falsely stated that the 
Horner method gives the smallest number of multiplications necessary for the evalu- 
ation of any polynomial (it might more appropriately have been described as usu- 
ally the shortest convenient method) ; in Chapter 7 she does not note that her matrix 
must have real eigenvalues, and so on. 

The book is not complete. The author does not suggest (for example) a double 
iteration for taking square roots (using the iteration noted above for the square 
root and a second iteration mentioned in her book for reciprocals so that no division 
is required). Only the power method of finding eigenvalues is mentioned, and, in 
general, only those methods for which codes are mentioned above are considered. 
There are frequent references to a companion volume on Automatic Digital Calcu- 
lators [2] and another on Numerical Methods [3], but these seem to the reviewer to 
be inadequate reference books upon which to base all the codes of a digital computer. 


C. B. T. 


1. C. Hastines, Approximations for Digital Computers, Princeton, 1955, [MTAC, v. 9, 1955, 
p. 121). 

2. A. D. Boorn & K. H. V. Boorn, Automatic Digital Calculators, 2nd edn., Butterworths, 
London, 1956. 

3. A. D. Bootu, Numerical Methods, 2nd edn., Butterworths, London, 1957, [MTAC, v. 10, 
1956, p. 105). 


30[X, Z).—Gusert R. Gray, DITMB Univac Transportation Simplex, Report 1266, 
David Taylor Model Basin, Washington, D. C., 1958, ii + 23 p. diagrs., tables, 
refs. 


This report describes the methods and UNIVAC I routines developed to obtain 
initial solutions for Suzuki’s Transportation Simplex Method. The initial solutions 
considered are the column, row, and matrix minimal (or maximal) and the North- 
west Passage solution. At present the m X n cost matrix is limited by m S 30 and 
m +n S 719. The starting routines developed serve two main purposes: (1) to 
obviate tedious hand calculations and data tape preparations, and (2) to reduce the 
machine time required to solve a complete transportation problem. 


AvuTHOR’s SUMMARY 


31{Z].—Rira M. Horsert, Artur SHapiro, Karen L. Brapuey, Sysit E. Jakos 
& Harriet R. Matin, Automatic Routines for Programming Management Data 
Problems on Univacs I and II, Report 1241, David Taylor Model Basin, Wash- 
ington, D. C., 1958, iv + 71 p. Tables. 


A collection of automatic programming routines is presented here. These routines 
are designed for the UNIVAC, a high-speed, automatic, electronic, digital computer. 
Routines are included which perform the functions of edit, sort, merge and data 
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conversion. This report contains a detailed description of the capabilities and limita- 
tions of the routines, and full directions on how to use them. 
AvuTHORS’ SUMMARY 
32(Z|.—Stpney Lecutrer, Survey of analog-to-digital converters, Report # 1257, 
David Taylor Model Basin, Washington, D. C., 1958, iv + 145 p. 


“Many test facilities are recording measurements of physical phenomena sensed 
by transducers having an analog output. The analog-to-digital converter can be of 
assistance to the data reduction field by providing the link, without human inter- 
vention, between these measurements and the modern high speed digital computer.” 
So states the author of this report. 

This report is a compilation of analog-to-digital converters, cataloged alpha- 
betically by manufacturers. It is to be used as a computer system engineers’ hand- 
book. The major contribution of this report is the readily available listing, grouping 
those analog-to-digital converters which meet the specifications of a particular 
design. Characteristics of each converter, or encoders as they are usually called, are 
listed with photographs. 

Listed characteristics are: power, space, weight, cost, input and output for- 
mats, and typical applications. 

In addition, a brief discription of the availability of these converters, their 
accuracies, additional features and remarks which may be pertinent to each encoder 
are included. A cross reference list of manufacturers and a list of conversion speeds 
are also incorporated in this report. 


ALEXANDER C. ROSENBERG 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, District of Columbia 


33[Z|—-MontTGomErRY Puister, Jr., Logical Design of Digital Computers, John 
Wiley and Sons, Inc., New York, 1958, vii + 408 p., 23 cm. Price $10.50. 


Logical design is a name given for the planning of the interconnection of cir- 
cuits, memories, and other units to form a digital system such as a computer. While 
a large amount of ingenuity and experimental thought may always be required for 
logical design, it is possible, nevertheless, for the designer to benefit greatly from the 
effective use of several mathematical techniques. Procedures utilizing some of these 
techniques are presented and well illustrated by examples and exercises. The pres- 
entation is complicated by the fact that these procedures must depend somewhat 
upon the nature of the circuits and other equipment being interconnected, but the 
author solves his problem by limiting his treatment, on the one hand, and by in- 
dicating procedural variations, on the other. Thus, the examples are principally of 
serial computers, and all are synchronous, or clocked. 

The author follows the methods of the “algebraic equation” school of logical 
design, although in his use of the Karnaugh-Veitch diagram throughout the process 
for representing switching functions he has made a significant change from the 
original approach. His discussion requires no specific background on the part of the 
reader, but a general familiarity with digital devices would be helpful as well as a 
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willingness to follow the many tedious details which are characteristic of logical 
design problems. The book should be useful not only to those who have had slight 
experience and who wish to learn logical design, but also to many who are accus- 
tomed to more haphazard and intuitive methods. 

In this “how to do it” book one finds that theory is given relatively light treat- 
ment; probably to conserve space and increase the number of potential readers. 
Therefore, the student wishing to deepen his understanding of the ideas which led 
to the development of the emphasized techniques may find it desirable to supple- 
ment his reading. The bibliography is also, understandably, slanted in the direction 
of practicality and one finds few references to the theory of automata as developed 
by Burks, Kleene, von Neumann, Turing, McCulloch, and others. True, such ma- 
terial would be of little value in the practical problems of logical design, but it is 
from this direction that the impetus for further mathematical development has 
been received. 


Davin E. MULLER 
University of Illinois 
Digital Computer Laboratory 
Urbana, Illinois 
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TABLE ERRATA 


269.—Paut F. Byrp & Morris D. FriepmMan, Handbook of Elliptic Integrals for 
Engineers and Physicists, Springer-Verlag, Berlin, 1954. 


I should like to point out an error appearing on page 340 of this book, where the 
value of KZ (8, k) corresponding to sin~' k = 75° and 6 = 20° is given erroneously 
as 0.565367 instead of the correct value, 0.565010. 


DomEnico CALico 
Instituto Nationale per le Applicazioni del Calcolo 
Rome 


Epiror1AL Nore: An independent calculation shows that the functional value 
in question should read 0.56501 09053 37977 to 15D. 
J. W. W. 


270.—Tosio Krracawa, Tables of Poisson Distribution, Baifukan, Tokyo, 1952. 


This table was originally reviewed in MT AC, April 1953, p. 92-93. At that time 
the reviewer found no discrepancies in about 100 comparisons with Molina’s table 
Subsequently in MTAC, April 1954, p. 95, D. Teichroew compared Kitagawa’s 
table with an unpublished table for m = 1(1)10. In the instance of Teichroew’s 
comparison, because of the large intervals between values of m checked, this 
amounted to comparison of 219 entries. 

The writers have made what is probably the most extensive comparison and 
check of this table to date. The work consisted of: 


1. For every m in the table, cumulating all entries to see if they added to 1 or 
closely so. This check will catch large errata; a number of ne. errata were 
so detected. 

2. For « = 0, for every m in the table, Kitagawa’s values were checked against 
NBS Applied Mathematics Series, No. 14, Tables of the Exponential Function 
e*, and the supplementary AMS-46, Table of the Descending Exponential. A 
few new rounding errors were thus disclosed. 

3. Roughly 300 entries were checked by direct computation, using the tables 
mentioned above and others. Over 200 of these were checked for various 
values of m < 2, using the table x"/n! in AMS-37, Tables of Functions and 
Zeros of Functions. One rounding error was detected and verified. 

For the convenience of the users of Kitagawa’s volume, the 5 errata noted by 


Teichroew have been incorporated in the present list. A number of poorly printed 
entries were also noted. 


m x Reads Should Read 

.305 6 .0000 0002 .0000 0082 
505 0 .-035 0558 .6035 0558 
.579 0 .5604 5855 .5604 5854 
.671 0 .5111 9713 .5111 9712 
.831 0 .4356 1346 .4356 1345 


141 





TABLE ERRATA 


3 
N 


Reads Should Read 
.0006 5050 .0006 5049 
.6000 0294 .0000 0294 
.0333 732- .0333 7327 
.1-67 7953 .1867 7953 
.0261 9892 .0361 9892 
.0223 479 .0228 479 
.1535 085 .1585 085 
.0024 083 .0025 083 
.0384 486 .0374 486 
.0902 262 .0912 262 
.1336 747 .1338 747 
.0649 185 .0659 185 
.0003 458 .0003 457 
.1369 0— 1369 060 
.0001 600 .0001 601 
.0194 306 .0194 307 
.0000 158 .0000 159 
10 .0204 075 .1204 075 
18 .0064 400 .0065 400 
11 .1137 363 1137 364 
15 .0347 180 .0347 181 
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Note: Hyphen indicates indistinct or missing figure. 


C. R. Sexton 
C. A. SEXTON 


J. A. SEXTON 
3009 Claremont Avenue 
Berkeley 5, California 


271.—Hans RIieEskEz, “Mersenne Numbers”, MTAC, v. 12, 1958, p. 207-213. 


On p. 210, the factor of 2? — 1 corresponding to p = 2689 should read 7158119 
instead of 7158199, and on p. 211 the factor corresponding to p = 5743 should read 
643217 instead of 543217. No other errors in this table were revealed by an in- 
dependent run on an IBM 704. 


J. L. SELFRIDGE 
International Business Machines Corporation 
Research Center 
Yorktown Heights, New York 





CORRIGENDA 


J. W. Wrencnu, Horace Scudder Uhler (1872-1956), MTAC, v. XII, 1958, 
p. 267: The title of item 28 in the bibliography should read, “‘Many-figure approx- 


imations for ¥/2, ¥/3, ¥/4 and ~¥/9 with x’ data.” 


J. W. W. 


K. A. Janrpov & 8. N. Razumovsk1, Tables of the logarithmic integral, MT AC, 
Review 106, v. XII, 1958, p. 238-239. 
for read 
p. 238, line 4 from the bottom, Izdatel’stro Izdatel’stvo 
p. 239, line 3, x S .5 zs 05 
p. 239, line 6, (1 — 7) (1 — ?t). 
A. F. 


E. N. Dexanosipze, Tables of cylinder functions of two variables, MTAC, 
Review 107, v. XII, 1958, p. 239-240. 
for read 
p. 240, line 13, decinal decimal. 
A. F. 


R. A. BuckincHaM, Numerical Methods, MTAC, Review 114, v. XII, 1958, 
p. 251. The publisher of this book has been listed as G. P. Putnam’s Sons. It should 
read Pitman Publishing Corporation. 











